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Abstract 

A  robust  and  efficient  optimization  code  is  developed  and 
validated.  The  code  is  used  to  redesign  an  existing  Mach  12  wind 
tunnel  nozzle  and  utilizes  response  surface  methodology  (RSM) 
techniques.  Explicit,  globally  second-order,  flux-difference¬ 
splitting  algorithms  are  used  to  solve  the  Navier-Stokes  (NS)  and 
Parabolized  Navier-Stokes  (PNS)  flow  solvers  incorporated  into  the 
optimizer  code.  Either  the  Baldwin-Lomax  or  the  Yang-Shih  k-e 
turbulence  model  may  be  employed  in  the  optimization  code. 

First,  2-D/axisymmetric  NS  and  PNS  flow  solvers  are 
developed/modified  and  account  for  perfect  gas/nonequilibrium 
chemically  reacting  flows.  All  solvers  are  validated  against 
Computational  Fluid  Dynamics  (CFD)  and  experimental  data. 

The  optimization  code  is  subsequently  developed  and  validated. 
The  optimization  code  is  then  used  to  optimize  the  Mach  12  nozzle 
design  and  the  computed  results  are  compared  with  those  of  the 
original  nozzle.  The  code  is  tested  for  robustness  and  on  three 
separate  occasions  locates  the  global  minimum  synonymous  with  the 
"global  best"  optimized  nozzle.  Though  an  optimized  nozzle  is 
obtained,  it  is  not  as  free  of  disturbances  in  the  uniform  inviscid 
core  at  the  exit  as  possibly  desired. 
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DEVELOPMENT  AND  TESTING  OF  A  NEW  OPTIMUM  DESIGN 


CODE  FOR  HYPERSONIC  WIND  TUNNEL  NOZZLES, 
INCLUDING  BOUNDARY  LAYER,  TURBULENCE, 

AND  REAL  GAS  EFFECTS 

1 .  Introduction 

This  chapter  presents  the  problem  addressed  by  the  current 
research,  a  short  history  of  previously  conducted  research,  the 
purpose  and  objectives  for  the  present  study,  and  an  overview  of 
the  approach  of  this  research  to  provide  a  foundation  for  the 
following  chapters  of  this  dissertation.  It  begins  with  a 
statement  of  the  problem  regarding  previous  and  current  design 
methods  of  hypersonic  wind  tunnel  nozzles,  which  have  inherent 
restrictions  on  accuracy  or  a  high  level  of  user  interaction 
requirements.  The  chapter  continues  with  a  presentation  of 
recently  conducted  research  in  the  area  of  hypersonic  nozzle  design 
and  research  on  the  testbed  nozzle  for  this  investigation;  the 
previously  accomplished  research  aids  in  validating  the  developed 
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code.  The  chapter  concludes  with  the  purpose  of  this  research  and 
concludes  with  an  overview  of  the  research  approach. 

1 . 1  Problem 

Both  previous  and  current  hypersonic  wind  tunnel  nozzle  design 
techniques  have  unique  issues  associated  with  them  that  result  in 
some  level  of  limitation  on  their  application.  The  method  of 
characteristics  and  boundary  layer  (MOC/BL)  design  technique  was 
used  by  the  hypersonic  nozzle  designers  of  yesteryear.  The  issue 
associated  with  this  technique  is  its  lack  of  accuracy  in  computing 
the  flow  in  a  hypersonic  nozzle.  Current  design  techniques  are 
more  accurate.  However,  in  such  methods,  the  optimization  scheme 
does  not  account  for  nonequilibrium  and  chemical  reaction  effects 
and  is  somewhat  restricted  in  its  application.  Such  methods  also 
require  a  high  level  of  interaction  and  knowledge  on  the  part  of 
the  design  engineer. 

In  spite  of  the  fact  that  hypersonic  wind  tunnels  have  been 
used  for  the  last  30  to  40  years  and  that  incredible  advances  have 
been  made  in  CFD  over  the  same  time  period,  most  hypersonic  wind 
tunnels  were  designed  based  upon  the  MOC/BL  method  developed  by 
Prandtl  and  Buseman  in  1929  [1] .  In  this  method,  the  inviscid  core 
is  determined  through  the  MOC  and  is  augmented  by  attaching  a 
displacement  thickness  obtained  from  a  BL  solution.  In  a  high  Mach 
number  wind  tunnel  nozzle  (M  >  8),  with  increasing  design  Mach 


2 


number  this  MOC/BL  method  suffers  an  increasing  loss  in  accuracy  as 
is  shown  below  [2] .  Such  a  loss  in  accuracy  may  not  be  acceptable 
for  today's  hypersonic  aerospace  vehicle  designers  or  rocket  nozzle 
designers;  wind  tunnel  nozzles  designed  based  on  the  MOC/BL  method 
may  be  delivering  test  conditions  different  than  those  predicted  by 
the  original  nozzle  solution.  Such  wind  tunnel  nozzles  are  not 
delivering  the  predicted  flow  to  the  wind  tunnel  test  section,  and 
therefore  designers  of  hypersonic  vehicles  are  using  less-than- 
desirable  test  section  flow  to  obtain  experimental  results  and  to 
aid  in  validating  their  CFD  solutions  with  experimental  results. 
The  hypersonics  community  has  recognized  this  shortcoming 
associated  with  existing  nozzles,  and  is  currently  focusing  much 
CFD  research  and  development  on  this  problem. 

Current  hypersonic  nozzle  design  procedures  are  much  more 
accurate  due  to  recent  advances  in  CFD,  but  still  fall  short  on 
several  major  fronts.  First,  the  current  optimization  codes  do  not 
account  for  full  vibrational  nonequilibrium,  chemical  reaction 
effects  which  would  likely  be  significant  in  a  hypersonic  nozzle. 
Second,  the  application  of  such  optimization  codes  requires  a 
smooth  expansion  corner  in  order  to  model  the  centerline  Mach 
number  distribution  which  is  used  to  drive  the  optimization  scheme. 
This  requirement  somewhat  limits  the  applications  to  which  the 
optimization  code  may  be  applied.  And  third,  successful 
interaction  at  the  necessary  level  with  such  codes  demands  that  the 
design  engineer  be  quite  knowledgeable  of  the  optimization 
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procedure;  due  to  the  complexity  of  the  optimization  scheme/code, 
the  required  understanding  of  the  optimization  procedure  is  at 
such  a  detailed  level  as  to  necessitate  intensive  training. 

There  is  therefore  a  current  need  within  the  hypersonics 
community  for  an  accurate,  hypersonic  nozzle  design  procedure  which 
accounts  for  the  possibility  of  chemical  reaction  effects,  is  broad 
in  its  potential  applications,  and  requires  a  minimum  of  computer 
time  and  designer  interaction. 

1.2  History  of  Nozzle  Research 

This  section  presents  two  main  topics:  a  history  of  research 
in  hypersonic  nozzle  design,  and  a  history  of  research  conducted  on 
the  testbed  for  this  investigation,  that  of  the  Wright  Laboratory 
(WL)  20"  Mach  12  wind  tunnel  nozzle. 

1.2.1  History  of  Hypersonic  Nozzle  Design.  This  section  is 
broken  up  into  two  parts  which  span  the  time  in  which  numerical 
methods  have  been  used  to  design  hypersonic  nozzles.  The  first 
part  covers  the  classical  MOC/BL  method  previously  mentioned.  The 
second  part  covers  more  recent  hypersonic  nozzle  design  methods. 

1.2. 1.1  The  MOC/BL  Method.  The  MOC/BL  method  is  an 
iterative  design  process  in  which  the  inviscid  core  obtained  by  the 
MOC  and  the  boundary  layer  interact  in  a  loosely  coupled  fashion 
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through  their  respective  boundary  conditions  (Fig.  1) .  Either  a 
wall  boundary  is  specified,  from  which  the  displacement  thickness 
and  inviscid  core  can  be  calculated,  or  the  inviscid  core  and 
displacement  thickness  are  calculated,  from  which  the  location  of 
the  wall  boundary  is  determined.  The  iterations  continue  until 
some  form  of  tolerance  has  been  reached  with  regard  to  the  location 
of  the  wall  or  some  other  parameter (s) . 

A  detailed  example  of  the  MOC/BL  process  is  given  in  detail  in 
[2]  .  In  this  particular  approach,  the  flow  inviscid  core  is  first 
generated  using  the  MOC.  The  displacement  thickness  is  then 
generated  using  the  BL  method  and  is  added  to  the  inviscid  core  to 
determine  where  the  "hard"  wall  is  located.  Boundary  layer 
corrections  are  subsequently  performed  and  the  resultant 
displacement  thickness  subtracted  from  the  "hard"  wall.  At  this 
point,  the  resulting  streamline  at  the  edge  of  the  displacement 
thickness  is  compared  to  the  inviscid  streamline  obtained  through 
the  MOC  (Fig.  1)  .  The  nozzle  exit  inviscid  core  diameter  is  then 
resized,  if  necessary,  to  obtain  a  better  match  between  the 
inviscid  core  and  displacement  thickness  edge  streamlines.  Note 
that  in  this  iterative  method  the  exit  diameter  is  constrained 
while  the  nozzle  length  is  not. 

As  previously  mentioned,  the  MOC/BL  method  becomes  inaccurate 
for  hypersonic  nozzles  with  a  design  Mach  number  of  8  or  more;  this 
inaccuracy  is  a  result  of  one  major  assumption  inherent  to  this 
method  [2]  :  the  region  of  flow  between  the  boundary  layer  edge  and 


5 


the  inviscid  streamline  is  assumed  to  be  inviscid  and  irrotational 


(Fig.  2)  .  However,  in  reality,  the  flow  in  this  region  is  both 
viscous  and  rotational,  causing  the  "characteristic  waves  to  curve 
toward  the  wall  as  the  Mach  number  decreases"  [2]  (Fig.  2)  .  Thus, 
the  waves  would  likely  not  be  reflected  as  they  are  assizmed  to  with 
the  inviscid,  irrotational  assumption.  This  deviation  is  small  for 
low  Mach  numbers  due  to  the  incidence  angle  of  the  characteristic 
waves,  which  results  in  less  distance  spent  in  the  outer  viscous 
portion  of  the  boundary  layer  (Fig.  3)  .  However,  for  high  Mach 
numbers,  the  incidence  angle  decreases,  causing  the  characteristic 
wave  to  "remain  in  the  viscous  region  for  a  much  longer  distance" 
[2]  .  Thus,  "the  effects  of  rotational  flow  would  be  more 
significant"  [2] .  Benton  et  al.  conclude  by  recommending  "that  a 
(hypersonic  nozzle)  design  method  using  the  Parabolized  Navier- 
Stokes  Equations  be  developed"  [2] . 

1.2. 1.2  History  of  Recent  Hypersonic  Nozzle 
Design/Optimization  Research.  Most  of  the  current  research  in  the 
optimization  area  has  come  about  through  the  work  of  Korte,  at  NASA 
Langley.  At  the  Naval  Surface  Warfare  Center's  Mach  14  Nozzle 
facility,  Korte  et  al.  have  studied  hypersonic  flow  physics  through 
both  computational  and  experimental  efforts  [9]  .  The  research 
focused  on  providing  better  understanding  of  the  flow  physics  for 
wind  tunnel  redesign  purposes.  The  computational  calculations 


6 


utilized  the  CAN-DO  code  developed  by  Korte  [9]  .  This  computer 
code  was  modified  to  account  for  vibrational  equilibrium  effects, 
and  used  the  Baldwin-Lomax  (B-L)  and  Cebeci-Smith  (C-S)  turbulence 
models.  The  vibrational  equilibrium  effects  assumed  thermally 
perfect,  calorically  imperfect  gas.  The  gas  that  was  studied  was 
nitrogen  (Nj) .  The  computational  results  agreed  well  with  the 
experimental  data. 

The  Korte  optimization  method  provided  good  results  for  the 
test  case  in  [9]  as  well  as  those  in  [60]  .  However,  as  previously 
mentioned  the  method  does  have  some  inherent  limitations.  In 
boosting  computational  efficiency,  the  code  was  designed  to  account 
for  at  most  vibrational  equilibrium  in  the  wind  tunnel  test  gas. 
The  possibility  to  design  a  hypersonic  nozzle  and  account  for 
vibrational  nonequilibrium  and  chemical  reaction  effects  is  not 
currently  available. 

Additionally,  the  CAN-DO  code  requires  a  macroscopically 
smooth  expansion  corner  in  order  to  utilize  a  modification  of 
Sivells'  method  for  modeling  the  centerline  Mach  number 
distribution  theoretically  (Figs.  4  and  5)  [9,36].  This 
distribution  is  used  in  the  optimization  scheme  to  drive  the  design 
nozzle's  computed  centerline  Mach  number  distribution  to  match  the 
theoretical  one  [9] .  As  previously  mentioned,  the  requirement  of 
a  smooth  expansion  corner  somewhat  limits  the  applications  to  which 
the  optimization  code  may  be  applied. 
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Finally,  the  CAN-DO  code  utilizes  a  complex  optimization 
scheme  which  requires  the  nozzle  wall  to  be  specified  with  cubic 
splines  sections.  The  nozzle  design  engineer  must  specify  the 
number  of  cubic  spline  knots  (junctions  or  nodes  between  adjacent 
cubic  spline  sections  which  define  the  nozzle  wall)  for  the  nozzle 
contour.  The  design  engineer  must  also  specify  the  axial  location 
of  the  control  points  associated  with  the  cubic  spline  knots  on  the 
nozzle  contour.  Determining  the  number  of  knots  on  the  cubic 
spline  and  the  location  of  the  control  points,  is  not  a 
straightforward  or  easy  process"  [9] .  Additional  knots  must  be 
"added  until  a  solution  is  obtained  which  satisfies  the  design 
requirements"  [9]  .  All  of  these  combine  into  a  requirement  for 
much  interaction  and  knowledge  on  the  part  of  the  design  engineer. 

In  further  research  associated  with  similar  methods,  Hackett 
has  computed  CFD  flow  field  solutions  for  the  NASA  Langley  Research 
Center's  15-Inch  Mach  6  High  Temperature  Tunnel  and  16-Inch  Mach  17 
Nitrogen  Tionnel  [26]  .  Full  NS,  PNS,  and  Method  of  Characteristics 
techniques  were  used,  and  the  solutions  were  compared  with  tunnel 
calibration  data.  The  research  was  performed  to  address  the  need 
for  high-quality  wind  tunnel  data  against  which  to  validate  CFD 
codes .  The  comparisons  between  the  CFD  solutions  and  the 
experimental  data  for  the  Mach  6  tunnel  were  shown  to  be  in  good 
agreement  (nozzle  exit  Pitot  pressures  up  to  ±3%  different  and 
excellent  agreement  in  the  exit  Mach  number  profiles  were  seen 
between  CFD  solutions  and  the  data) .  However,  the  comparisons  for 
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the  Mach  17  tunnel  were  not  (nozzle  exit  Pitot  pressures  showed 
differences  of  up  to  20%  difference  between  CFD  solutions  and 
experimental  data) .  This  was  attributed  largely  to  inaccuracies  of 
turbulence  models  for  hypersonic  flow  (the  B-L  turbulence  model  was 
based  upon  incompressible  low  speed  flows) .  Degradation  of 
accuracies  at  higher  Mach  numbers  was  expounded  upon,  and  the 
possibility  of  the  necessity  of  compressibility  effects  and 
corrections  was  proposed.  Also  identified  as  a  factor  for  the  poor 
agreement  was  the  fact  that  nitrogen  is  not  calorically  nor 
thermally  perfect  at  such  high  temperatures  and  pressures. 

Several  research  efforts  have  concentrated  on  the  need  to 
include  non-ideal  gas  effects  in  the  full  NS  solution  of  high  Mach 
nozzle  flow.  For  instance,  in  the  work  done  by  Johnson  et  al.,  it 
was  shown  (using  the  MOC/BL  method)  that  an  increase  in  the 
displacement  thickness  of  13%  occurred  when  high  temperature  gas 
effects  were  included  as  compared  with  perfect  gas  results  [30] . 
In  addition,  the  inviscid  core  exit  radius  increased  9%  when  high 
temperature  effects  were  included.  It  should  be  noted  that  the 
stagnation  conditions  for  this  research  were  somewhat  extreme 
(stagnation  pressure  of  1.01325x10®  N/m^  and  stagnation  temperature 
of  2333°K)  [30] .  However,  the  data  appears  to  support  the  need  to 
investigate  whether  high  temperature  gas  effects  will  have  a 
significant  impact  on  a  particular  high  Mach  nozzle  flow. 

In  a  high  Mach  nozzle  study  performed  by  Candler  and  Perkins, 
a  number  of  interesting  facts  relating  to  non-ideal  gas  effects 
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were  presented  [31] .  First,  for  the  hypersonic  nozzles  tested, 
vibrational  nonequilibrium  provided  its  major  effect  upstream  of 
the  throat;  the  vibrational  temperature  froze  near  the  throat  and, 
downstream  of  the  throat,  the  flow  essentially  acted  as  a  perfect 
gas  with  constant  y=1-4.  Second,  vibrational  nonequilibrium  and 
perfect  gas  solutions  gave  very  similar  results,  while  equilibrium 
flow  provided  a  substantially  lower  Mach  number  at  the  nozzle  exit. 
And,  third,  the  zero-equation  algebraic  turbulence  model  used  has 
a  significant  effect  upon  the  flow  solution  obtained  (C-S,  B-L,  and 
Renormalization  Group  (RNG)  models  were  used)  ;  a  2%  to  4% 
difference  was  seen  in  the  inviscid  core  exit  Mach  number  when  the 
different  turbulence  models  were  used. 

To  summarize,  the  history  of  recent  hypersonic  nozzle  design 
and  optimization  points  to  a  number  of  critical  hypersonic  nozzle 
design  community  needs.  One,  an  accurate,  efficient,  simpler 
optimization  code  would  be  useful.  Two,  a  given  nozzle  should  be 
testable  as  a  minimum  with  a  nonequilibrium  CFD  code  to  determine 
if  nonequilibrium  effects  need  to  be  accounted  for  in  the 
optimization  procedure.  If  nonequilibrium  effects  prove  to  be 
significant,  a  nonequilibrium  code  should  be  an  option  such  that  it 
can  at  the  very  least  be  used  in  the  final  iterations  of  the 
optimization  procedure.  Three,  the  turbulence  model  used  has  a 
significant  effect  on  the  nozzle  solution  and  its  selection, 
therefore,  must  be  addressed  in  a  critical  manner.  These  three 
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required  critical  technologies  formed  the  basis  of  the  present 
research. 

1.2.2  The  Wright  Laboratory  Mach  12  Nozzle  -  The  Test  Case. 

Since  data  and  personal  expertise  were  locally  available  for  the 
Wright  Laboratory  (WL)  20"  Mach  12  wind  tunnel  nozzle,  this 
facility  was  chosen  to  be  the  testbed  for  this  investigation.  The 
Air  Force  WL  hypersonic  wind  tunnel  is  an  intermittent,  blowdown- 
to-vacuum  facility  with  an  axisymmetric  20  inch  exit  diameter 
nozzle  contoured  to  produce  uniform  Mach  12  or  14  flow  in  an  open 
jet  test  section  (Fig.  6)  .  For  the  Mach  12  nozzle,  the  tunnel 
reservoir  pressure  ranges  from  4.3169x10®  to  11.0316x10®  N/m^  (600 
to  1600  psia) ,  the  Reynolds  number  ranges  from  1.3123  to  3.2808 
million  per  meter  (0.4  to  1.0  million  per  foot),  the  stagnation 
temperatures  range  from  1000  to  1111  degrees  Kelvin  (1800  to  2000 
degrees  Rankine) ,  and  run  times  vary  from  5  to  8  minutes  [63]. 

A  number  of  studies  have  been  performed  on  the  WL  Mach  12 
nozzle.  Buck  et  al.  performed  experiments  which  showed  the  flow  to 
be  isentropic,  at  least  in  the  high  Mach  core  [3]  .  Trollier  et 
al.,  of  Science  Applications  International  Corporation  (SAIC)  under 
contract  with  the  Air  Force  Wright  Aeronautical  Laboratories,  used 
CFD  techniques  based  on  the  full  NS  and  PNS  equations  to  solve  the 
flow  fields  of  the  wind  tunnel  throat  and  supersonic  nozzle 
respectively.  It  was  demonstrated  computationally  that  the 
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experimental  exit  Pitot  profile  and  the  exit  Mach  profile  were 
bracketed  by  the  fully  laminar  CFD  profiles  on  the  one  side  and 
fully  turbulent  CFD  profiles  on  the  other  (Fig.  7)  [4] .  It  is  this 
author's  opinion  that  this  lack  of  agreement  with  the  measured  exit 
data  was  likely  due  to  uncertainty  in  the  transition  location, 
inaccuracies  in  the  turbulence  modelling,  or  simplified  modeling  of 
the  physics  of  the  flow  (that  is,  not  accounting  for  nonequilibrium 
effects)  . 

As  to  the  use  of  this  facility  as  the  testbed  in  the  present 
research,  it  is  the  author's  intention  to  use  both  experimental 
data  from  and  CFD  solutions  to  the  WL  Mach  12  nozzle  to  validate 
CFD  solutions  obtained  herein.  Additionally,  experimental  data 
from  the  WL  Mach  12  nozzle  will  be  used  to  compare  with  the  newly 
designed  Mach  12  nozzle  obtained  with  the  optimization  procedure. 

1.3  Purpose  and  Objectives  of  Present  Research 

The  purpose  of  this  research  is  to  redesign  an  existing 
hypersonic  wind  tunnel  nozzle  and  in  the  process  develop  a  new 
design  tool  with  the  necessary  inherent  methodology.  In  the 
process  of  redesigning  the  nozzle  it  will  be  desired  that  the  new 
design  provide  the  largest  uniform  flow  test  section  possible  with 
a  minimum  of  computer  time  and  designer  interaction.  To  achieve 
this  end,  state-of-the-art  turbulence  modelling  and  axisymmetric 
capabilities,  viscous  effects,  and  compressibility  effects  must  be 
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accounted  for  in  both  NS  and  PNS  codes. 


In  addition,  an 


optimization  scheme/code  is  required  to  perform  nozzle  design.  The 
computer  program  developed  will  be  an  improved  tool  for  the 
redesign  of  existing  MOC/BL-designed  wind  tunnels  due  to  the 
complexity  of  the  flow  solvers  and  the  simplicity  of  the 
optimization  scheme.  The  developed  optimization  code  can  later  be 
used  for  the  design/redesign  of  more  accurate  hypersonic  wind 
tunnel  nozzles  or  other  high  Mach  propulsive  devices. 

Thus,  there  are  three  objectives  to  this  research,  one  primary 
and  two  minor  but  valuable  secondary  objectives: 

-  The  primary  objective  is  to  develop  and  validate  a  computer  code 
to  be  used  in  the  optimal  redesign  of  the  testbed  nozzle.  The 
code,  with  state-of-the-art  NS  and  PNS  solvers  incorporated,  will 
accurately  and  robustly  select  an  optimal  set  of  design  parameters 
based  on  some  established  design  criteria,  utilizing  a  simple 
optimization  scheme.  Specifically,  the  code  will  be  used  to  design 
a  hypersonic  nozzle  with  the  constraints  of  fixed  length,  throat 
radius,  and  exit  radius  based  on  the  main  application,  the  test 
case  described  in  Sect.  1.2.2. 

-  A  secondary  objective  is  to  determine  if  the  B-L  turbulence  model 
or  the  Yang-Shih  (Y-S)  k-e  model  would  help  provide  a  more  accurate 
solution  to  the  flow  field  of  the  testbed  nozzle. 

-  The  other  secondary  objective  is  to  test  the  newly  optimized 
nozzle  for  off-design  conditions.  Specifically,  the  new  nozzle 
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design  will  be  tested  with  changes  in  back  pressure  and  the  results 
analyzed. 

1.4  Overview  of  Dissertation  Approach 

In  order  to  achieve  the  objectives  mentioned  above  a  number  of 
critical  steps  needed  to  be  taken.  The  first  step  taken  was  to 
develop  and  validate  a  perfect  gas  NS  axisymmetric  laminar  code. 
The  second  step  was  to  develop  and  validate  a  perfect  gas  and  a 
nonequilibrium,  chemically  reacting  (NECR)  NS  axisymmetric 
turbulent  code.  With  the  preceding  steps  completed,  a  comparison 
could  be  made  between  nozzle  flow  solutions  obtained  with  the  two 
turbulence  models  incorporated  in  the  second  step;  thus  a 
turbulence  model  could  be  selected  for  the  remainder  of  the 
research,  satisfying  one  of  the  secondary  objectives.  The  third 
step  entailed  developing  and  validating  a  perfect  gas  PNS 
axisymmetric  laminar/turbulent  code.  The  fourth  step  was 
development  and  validation  of  a  simpler  optimization  scheme  and 
code,  with  siabsequent  redesign  of  the  testbed  nozzle,  and  tests  for 
robustness  performed  on  the  optimization  code.  Having  completed 
all  of  the  preceding  steps  would  satisfy  the  primary  objective  of 
this  research.  The  fifth  and  final  step  was  to  run  the  off-design 
cases  of  the  newly  designed  nozzle  to  determine  the  effects  of  a 
difference  in  back  pressure.  This  last  step  would  satisfy  the 
other  secondary  objective. 
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2.  Governing  Equations  and  Supporting  Theory 


The  equations  of  fluid  mechanics  as  pertain  to  the  various 
developed  codes  are  presented  in  this  chapter  along  with  important 
principles  of  optimization.  First,  the  general  form  of  the  NS 
axisymmetric  equations  is  given  in  Cartesian  and  computational 
coordinates,  followed  by  the  perfect  gas  and  NECR  forms.  Next,  the 
perfect  gas  PNS  axisymmetric  equations  are  presented.  Then  the 
turbulence  model  formulas  are  given.  Finally,  the  pertinent 
methodology  used  in  the  optimization  scheme  is  presented. 


2.1  NS  Axisymmetric  Equations,  General  Form,  Cartesian  Coordinates 

As  mentioned  previously,  the  flow  field  of  interest  is  that  of 
a  hypersonic  axisymmetric  wind  tunnel  nozzle.  Thus,  the  axi¬ 
symmetric  form  of  the  NS  equations  are  used.  The  equations  in 
vector  form  for  axisymmetric  coordinates  are  [48] : 


U*F  *  G  -  b 

t  X  y 


{H-G) 


in  non-conservation  form,  and 


-(yG)  -  5- 

^  y  J  y 


in  strong  conservation  form,  where 


U  -  (p,  pu,  pv,  E)' 


(1) 


(2) 


(3) 
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(4) 


F-  (pu,  puv-T^y,  u(E*p) 


G-  (pv,  puv-T^^,  pv^p-T^^,  viE*p)-uX^^-VT^^*q^)^  (5) 

H  =  (0,  0,  P-Tgg,  0)^  (6) 

(7) 

T^y  =  (8) 

Tyy  =  (2v^-u^-5^)  (9) 

Tee  =  (-i^;c-V2^)  (10) 

-{k^*k^)T^  (11) 

q^-  .(k^*kt)T^  (12) 


(Note  that  Stokes'  hypothesis  has  already  been  invoked;  ie, 
A=-2/3p) . 

2.2  NS  Axlsymmetrlc  Equations,  General  Form,  Computational 
Coordinates 

The  form  of  the  equations  given  above  is  for  that  of  a 
Cartesian  coordinate  system.  A  generalized  transformation  from  the 
physical  to  the  computational  domain  is  made  ( (x, y) <=> ( ^ , p ) )  in 
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order  to  facilitate  the  CFD  computations,  where  ^  =  ^{x,y)  is  the 
streamwise  computational  coordinate,  and  r|  =  ri(x,y)  is  the 
crossflow  computational  coordinate.  In  the  new  coordinate  system, 
Eq.  2  now  becomes  [48] 


U^*F,*d  -  6  — 
t  5  n  j 

where 


(13) 


U  -  yU/J 

(14) 

F= 

(15) 

G  -  ^(n^F^ri^G) 

(16) 

and  J,  the  coordinate  transformation  Jacobian,  is  defined  as 


(17) 


2.3  Extension  of  NS  Axisyrametric  Equations  to  Perfect  Gas  Form 

The  extension  of  the  governing  equations  from  general  form  to 
the  perfect  gas  form  are  now  shown.  First,  the  pressure,  p,  in  a 
thermally  and  calorically  perfect  gas  is  given  by: 

p  =  (Y-l)pej  (18) 
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while  the  total  energy  is: 


E  -  p(e_j+ 


2 


v")) 


(19) 


where  e^  =  c^T  (c^  a  constant)  for  a  thermally  and  calorically 
perfect  gas. 

The  molecular  viscosity  is  obtained  with  Keye's  Law 


j.l/2 


where 


::  1.490x10'®  N.sec/m^ 


(20) 


(21) 


=  122.22°K 


(22) 


and 


=  5°K. 


(23) 


The  molecular  thermal  conductivity  can  be  obtained  through  the 
use  of  the  Prandtl  number,  with 


(24) 


and  the  turbulent  thermal  conductivity  is  obtained  through  the  use 
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of  Pr^,  through 


k 


t 


Jiifp 

Pr^ 


(25) 


once  the  eddy  viscosity  is  known  from  the  turbulence  model. 

For  a  perfect  gas,  the  speed  of  sound,  a,  is  defined  as 

-3^  -  YP/P  -  (26) 

which  is  the  same  as  the  frozen  speed  of  sound,  since  in  both  cases 
it  is  assumed  that  no  chemical  reactions  are  occurring.  However, 
in  the  perfect  gas  case,  y  is  assumed  constant;  in  the  frozen  case, 
Y  may  be  a  function  of  temperature. 


2.4  NS  AxisynaaGtrlc  Equations ^  NECR  Form 

The  governing  equations  presented  to  this  point  have  dealt 
with  the  perfect  gas.  case,  but  to  determine  if  it  were  necessary 
to  account  for  high  temperature  effects,  it  was  necessary  to  use 
and  modify  the  computer  code  of  Josyula  [40]  .  This  code  accounts 
for  chemical  reactions,  as  well  as  nonequilibrium  vibrational 
effects.  The  form  of  the  governing  equations  is  supplemented  to 
account  for  five  individual  species  (O2,  Nj,  N,  0,  NO)  and ’three 
additional  energy  equations  which  track  the  vibrational  energy  of 
the  three  species  of  diatomic  molecules  (O2,  N2,  NO)  ,  for  a  total 
of  eleven  conserved  variables  at  any  given  point. 
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To  illustrate,  the  U  vector  from  Eq.  2  becomes 


Po/  Po'  Pm'  P^'  P^'  Pj.,%'  Po,%'  Pm^y„/ 


(27) 


where  p,  the  total  density,  is  now  equal  to  the  sum  of  the 
individual  species'  densities  -  Po2*Ph^*Po*Pn*Pw^  •  ^  vector 

from  Eq.  2  is  now 


F  - 


V 

Po,“ 

P»r“ 

Po“ 

Pj^o“ 

pu^*p-X^ 

e  u+g 

^0^  Vr 


-2  '0^ 

Pm^v 

^HO  ^NO 


(28) 


the  G  vector  from  Eq.  2  is 


G  - 


P,r3^ 

Po/ 

Pn^ 

Po^ 

Pm^ 

puv-T 


yx 


pv'^+p-T 


yy 


NO 


NO 


{E*p-T  )  V-T  u*q 


(29) 
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and  the  H  vector  from  Eq.  2  becomes 


where  the  total  energy,  E,  from  Eq.  19  has  become 


(30) 


E  -  De„  —  D(u^  +  v^)+D.,  e  ♦D«e  +D.,^e 

^  Tr-fRot  2^'  '  ^■*^2  ’'O  ’^HO 


The  additional  terms  to  the  above  vectors  account  for  the  chemical 
reactions  which  are  occurring  and  for  energy  being  transferred  from 
the  translational  energy  mode  to  vibrational  energy  as  a  result  of 
molecular  collisions  [50] . 

The  chemical  reaction  source  terms  for  the  ith  species,  in 
equation  29  derive  from  the  law  of  mass  action  and  consider  the 
five  significant  reactions  which  may  occur  between  the  different 
species  of  the  air  mixture  [40] .  These  possible  reactions  are 
represented  by 
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(32) 


^2 

M 

2N  *  M 

+ 

M 

20  *  M 

NO 

+ 

M 

f* 

N  *  0  *  M 

^2 

0 

NO  *  N 

NO 

0 

f* 

o 

to 

* 

in  which  M  is  merely  a  collision  body  such  as  an  atom  or  molecule. 

The  equilibrium  vibrational  energy  for  each  constituent  is 
modelled  using  the  simple  harmonic  oscillator  model  represented  by 
[33] 


e  - 


exp 

T 

-1 

1 

species. 

T 

V 

(33) 


where,  for  the  ith  species,  is  the  vibrational  temperature 


and 


0^^  is  the  characteristic  temperature  of  vibration  and  is  defined 


as  [33] 


0  z  JlLi 

k 

The  individual  factors  in  the  last  equation  are  Planck's  constant 
(h) ,  Boltzmann's  constant  (k),  and  the  frequency  of  oscillation  of 
the  ith  species  (v^)  . 

In  order  to  correctly  account  for  creation  of  vibrational 
energy,  the  source  terms  have  been  added  to  the  H  vector.  These 
source  terms  are  modelled  using  the  Landau-Teller  model  for  local 
relaxation  time  and  are  defined  as  [33] 
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e'*'  ( r)  -e 

''1 


(35) 


e 


T 


where  (r)is  the  vibrational  energy  associated  with  temperature  t 
for  the  ith  species,  and  is  the  local  relaxation  time  of  the 
ith  species  and  is  a  function  of  temperature  and  pressure.  The 
local  relaxation  time  for  O2,  N2  ,  and  NO  are  obtained  individually 
with 


T 


(36) 


where  x,  is  the  mole  fraction  of  the  ith  species  and  is  the 

interspecies  relaxation  time  between  species  i  and  species  k  and  is 
computed  with  the  form  developed  by  Millikan  and  White  [33]  . 
Vibration-vibration  coupling  is  assumed  negligible  and  not 
accounted  for.  However,  vibration-dissociation  for  the  diatomic 
species  is  computed  with  the  two-temperature  model  used  by  Park 
[51]  : 

-  ro-’ro-3  (37) 

in  which  is  the  dissociation  temperature  and  is  the  vibration 
temperature . 
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The  additions  to  the  viscous  portion  of  the  F  and  G  vectors 
are  the  vibrational  heat  fluxes,  which  are  given  by  [31,32] 


q 


^1 


(38) 


Even  though  the  perfect  gas  equation  of  state  no  longer  holds, 
the  new  equation  of  state  is  still  given  by  [35] 

p  ::  pRT  ,  (39) 


where 


(40) 


c 


p 


(41) 


c 


V 


(42) 


(43) 


Note,  however,  that  the  specific  heats  correspond  to  the 
equilibrium  portion  of  the  internal  energy,  i.e.,  the  translational 
and  rotational  energy. 

Finally,  the  last  equation  needed  is  that  for  the  speed  of 
sound.  Conveniently,  the  form  doesn't  change  from  that  of  the 
perfect  gas.  In  other  words,  the  equation  remains  [36] 
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(44) 


a  r  [  YP]1/2  ::  [yiTr] 

P 

As  pointed  out  in  [36],  "this  result  is  not  approximate,  but 
(actually)  corresponds  to  the  frozen  speed  of  sound  for  this  . . . 
nonequilibrium  flow." 


2.5  PNS  Equations,  Perfect  Gas  Form 

The  PNS  equations  are  a  modification  of  the  NS  equations,  but 
they  are  restricted  in  use  due  to  two  requirements:  "the  inviscid 
outer  region  of  the  flow  (must)  be  supersonic  and  the  streamwise 
velocity  component  (must)  be  everywhere  positive"  [48]  .  However, 
due  to  the  transformation  from  NS  to  PNS,  the  PNS  space-marching 
technique  allows  for  a  solution  to  be  obtained  in  the  same  time  as 
on  the  order  of  an  Euler  solution  [48]  .  As  shown  below,  space¬ 
marching  comes  about  as  a  result  of  dropping  the  time  derivative 
and  some  of  the  stress  terms  in  Eq.  13  and  integrating  in  the  ^  (or 
streamwise)  direction.  This  efficient  integration  scheme  is  what 
makes  the  PNS  equations  so  favorable  in  a  well-behaved  non- 
streamwise  separated  flow  [48] . 

Similar  to  the  NS  equations,  but  with  the  time  dependency 
removed  and  some  modifications  to  the  flux  vectors,  the  primary 
equation  in  transformed  coordinates  ,  similar  to  Eq.  13,  is  [49] 
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fO 

*F 

*G 

5 

1 

The  terms  in  the  brackets  on  the  left  hand  side  of  Eq.  45, 

are  the  Geometric  Conservation  Law  (GCL) 
terms  which  are  necessary  in  this  formulation  due  to  the  potential 
variation  of  the  grid  as  the  solution  marches  downstream  [9,49] . 
Note  that  the  ^  variable  is  only  a  function  of  x  (^  =  ^(x))  and  not 
of  y;  this  will  aid  in  the  transformation  from  the  computed 
variables  {F*)  to  the  conserved  variables  (U)  ,  as  will  be  shown 
below.  Note  also  that  the  Q  vector. 


Q  =  (0,  CO,  0,  0) 


(46) 


which  Korte  recommends  using  has  been  dropped  due  to  its  de¬ 
stabilizing  effect  on  the  PNS  solver,  as  recommended  by  White  et 
al.  [37]. 

In  the  flux  vectors,  F  and  G,  the  streamwise  (O  derivatives 
are  dropped,  whereas  in  the  F*  vector,  the  stress  terms  are 
dropped  completely.  Also,  in  F"^,  the  pressure  term  has  been 
changed  to  account  for  upstream  pressure  in  the  boundary  layer. 
This  change  comes  about  by  multiplying  the  pressure  by  Vigneron's 
coefficient  [49],  producing 

F*  -  (pu,  puu+cop,  puv,  (F*p)u)^.  (47) 


In  order  to  calculate  Vigneron's  coefficient,  co,  the  following 
ecjuation 
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(48) 


CO  = 


1* 


is  first  computed,  where  the  streamwise  Mach  number  is  denoted  by 
M^.  Then,  Vigneron's  coefficient  is  simply 

CO  -  min (1, ao)  (49) 

where  a  is  a  safety  factor  used  to  ensure  the  eigenvalues  of  the 
inviscid  portion  of  the  PNS  ecjuations  have  real  values  and  that  the 
equations  remain  hyperbolic  [49] .  The  value  used  for  a  in  this 
application  was  0,75.  This  value  seemed  to  provide  good  stability 
with  the  scheme  used.  Due  to  the  necessity  of  the  safety  factor, 
it  was  necessary  to  ensure  that  the  Mach  number  in  the  inviscid 
core  was  above  M ,  ,  which  is  defined  as 

>  [1/(ctY-Y*1)]1/2.  (50) 

With  a  having  the  value  of  0.75  and  assuming  a  Y  of  1*4,  the  value 
for  was  approximately  1.24. 

2. 6  The  Baldvin-Lomax  Turbulence  Model  Formulation 

Of  the  two  primary  zero-ecjuation  models  (B-L  and  C-S)  [28,  62], 
the  B-L  model  is  preferred  by  this  author  for  two  reasons.  First, 
one  need  not  locate  the  edge  of  the  boundary  layer  in  the  B-L 
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model.  Second,  the  B-L  has  been  validated  against  such  flows  as 
boundary  layer  over  a  flat  plate  and  transonic  airfoils  with 
excellent  comparative  results. 

Being  a  zero-equation  model,  the  B-L  turbulence  model  is  based 
upon  several  concepts  and  assumptions.  First,  the  necessary  length 
and  velocity  scales  are  obtained  from  the  mean  flow;  that  is, 
closure  for  the  eddy  viscosity  is  obtained  from  the  mean  velocity 
field.  Second,  it  is  assiamed  that  an  equilibrium  situation  exists 
between  the  mean  flow  and  the  turbulence  present.  Third,  for  the 
inner  region  of  the  boundary  layer  it  utilizes  a  form  of  Prandtl's 
Mixing  length  foinnulation,  which  defines  the  eddy  viscosity, 


dy 


(51) 


Fourth,  the  form  utilized  for  the  outer  region  of  the  boundary 
layer  is  modelled  after  the  Clauser  formulation  [28] . 

This  model  calculates  the  eddy  viscosity  (Pt)  by  calculating 
a  value  for  the  area  near  the  wall,  giving  (Pt)innerf  and  another 
value  for  the  flow  beyond  the  inner  region  (Pt)  outers  much  like  the 
classical  model  of  the  boundary  layer.  In  the  inner  region,  the 
eddy  viscosity  is  calculated  as 


=  pk  ^  -eXp  ( -^  )  ]  ^  |cd| 
A 

where 


(52) 
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CO  = 


du  dv 
dy  dx 


(53) 


(54) 


u 


T 


(55) 


with  k  and  constants  depending  on  the  flow  (originally  assigned 
values  of  0.0168  and  26,  respectively)  [28]. 

For  the  outer  region,  the  eddy  viscosity  is  defined  as  [28] 

outer  ~  ^  0.0168p)(1.6F^^,^) 

^kleb  (56) 


where 

F  .  =min(y  F  ,  0.25y  u.,.}/F  ) 

wake  '•'max  max'  •'max  diff  max' 

and 

^max  “  max(F(y) )  -  max(y|co|  [l-exp(-— )  ] ) 

A* 


(57) 


(58) 


where  y^^^  is  that  value  of  y  (for  a  given  x  value  held  constant) 
where  F(y)  =  F^^^  and 

^diff  ~ 
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=  [1*5.5(0.3-^)]-i  (60) 

•^max 

The  point  where  the  transition  from  the  inner  region  to  the  outer 
region  occurs  is  that  point  (in  the  y-direction  for  a  given  x- 
coordinate)  where  occurs  (y„ax)  • 

Generally  speaking,  the  B-L  model  produces  excellent  results 
in  flows  in  which  there  is  no  adverse  pressure  gradient.  However, 
it  is  well  known  that  it  doesn't  produce  a  reliable  solution  in 
separation  regions  and  recirculation  regions.  Hence,  it  is 
necessary  to  be  prudent  when  implementing  an  algorithm  for  this 
turbulence  model,  particularly  since  as  the  back  pressure  is 
increased  the  flow  might  develop  an  adverse  pressure  gradient  and 
separate . 

2.  7  The  Yang-Shih  k-e  Turbulence  Model  Formulation 

Originally  developed  by  Jones  and  Launder  [17]  to  compute 
incompressible  boundary- layer  flows,  the  two-equation  k-e  model  has 
become  a  reliable,  often-used  turbulence  model  in  the  CFD 
community.  A  comparison  of  numerous  2-equation  models  was  made  by 
Lang  and  Shih  [18] .  The  cases  tested  were  low  Reynolds  number,  2-D 
flows;  one  case  was  the  flow  over  a  flat  plate  and  the  other  was 
the  flow  through  a  fully  developed  channel.  The  Y-S  k-e  model 
appeared  to  be  one  of  the  most  accurate  and  robust  turbulent 
models,  and  it  was  one  of  the  few  models  which  approached  the 
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standard  k-e  model  away  from  the  wall.  For  these  reasons  this 
model  was  chosen  as  the  particular  two-equation  model  to 
incorporate  into  the  computer  program  developed  in  this  research. 

In  addition  to  the  standard  Y-S  k-e  model  as  proposed  by  Yang 
and  Shih  [52],  it  is  necessary  to  add  Sarkar's  compressibility 
correction  [19] .  This  necessity  comes  about  since  a  good  portion 
of  the  flow  in  a  hypersonic  nozzle  will  be  in  the  high  Mach  number 
region,  where  compressibility  effects  are  significant. 

There  are  two  additional  equations  which  must  be  calculated 
throughout  the  flow,  for,  this  being  a  two-equation  model,  there 
are  two  conserved  variables  needed  to  compute  the  eddy  viscosity. 
One  of  these  equations  governs  the  conservation  of  the  turbulent 
kinetic  energy,  k  (specifically  pk  is  conserved) ,  and  the  other 
governs  the  conservation  of  the  turbulent  energy  dissipation  rate, 
e  (again,  specifically  pe  is  conserved) . 

The  turbulent  kinetic  energy  conservation  equation  is  defined 
as 


dpk 

dt 


[puk-(ii+vi£/'7jt) 
ox  '  *  3x 


^  [pvk-(}a*p  /a  )  1^] 
ay  ay 


(61) 


and  the  turbulent  energy  dissipation  rate  conservation  equation  is 
defined  as 
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(62) 


9pe 

dt 


where 
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\  36  ) 

(72) 

exp(-0.004yU5e-^(y^)^2e-®(y^>^-8e-®{y^)^) 

(73) 

Re,  . 

(74) 

ve  jie 

-  0.09 

(75) 

Cel  =  1-44 

(76) 

C,,  =  1.92 

(77) 

a,  -  1.0 

(78) 

Ce  =  1*3 

(79) 

Pr,  =  0 . 90 

(80) 

The  compressibility  correction  already  appears  in  Eq.  61  as  the 
term  peoMt^,  and  represents  the  additional  dissipation  of  turbulent 
kinetic  energy  due  to  dilatation  (compressibility) ,  which  only 
appears  to  be  needed  for  hypersonic  flows  [20] . 

Note  that  an  additional  term  must  be  added  to  Eqs.  7  and  9  as 
a  result  of  the  dependence  of  the  eddy  viscosity  on  the  turbulent 
kinetic  energy,  k.  This  is  due  to  the  fact  that  the  Reynold's 
shear  stress,  -<UiUj>,  is  now  defined  as 
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whereas  before,  the  Reynold's  stress  had  only  been  modeled  by  the 
first  term  on  the  right  hand  side  of  Eq.  81.  The  terms  in  Eq. 
81  represent  the  mean  velocity  components  of  the  flow,  while  the 
terms  represent  the  fluctuating  velocity  components.  In 
particular,  -2/3pk  must  be  added  to  each  of  Eqs.  7  and  9,  producing 

and 

Tyy  =  (83) 

For  the  perfect  gas  case,  assuming  that  the  NS  solver  were  to 
be  used,  Eq.  2  would  look  the  same  as  originally  presented  but  with 
one  additional  source  term  on  the  right  hand  side,  producing 

-5-*S  (84) 

^  ^  y  J  y 


where  S  is 

5-  (0,  0,  0,  0,  P^-pe{l*(m^^)*D,  (85) 

■‘t 


The  dependent  variable  vector,  U  and  the  flux  vectors,  F  and  G,  now 
include  the  calculations  for  pk  and  pe.  To  illustrate,  the 
following  changes  are  made:  Eq.  3  becomes 


(86) 


U-  (p,  pu,  pv,  E,  pk,  pe)^, 

with  Eqs .  4  and  5  becoming 
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Additionally,  the  total  energy  now  accounts  for 
kinetic  energy,  as  well  as  turbulent  kinetic  energy, 
represented  by 


(87) 


(88) 


internal  and 
and  is 
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(89) 


B  =  p  (e^4.-|  {u^*v^)  *k)  . 

Inclusion  of  the  turbulent  kinetic  energy  into  the  total  energy  was 
done  per  the  work  of  Morrison  [39]  and  has  been  shown  to  aid  in 
convergence  of  the  solution. 

2.8  Supporting  Theory  for  the  Optimization  Scheme 

A  new  wind  tunnel  design  is  obtained  efficiently  by 
formulating  an  optimization  problem.  This  must  be  done  such  that 
the  end  result  is  a  wind  tunnel  nozzle  which  has  maximized  the  area 
of  the  inviscid,  high  Mach  core  at  the  exit  plane  with  minimal 
crossflow  and  as  few  as  possible  disturbances  in  that  core.  In 
other  words,  the  solution  to  the  optimization  problem  results  in  a 
smooth  expansion  and  uniform  flow  at  the  exit  plane. 

Three  different  optimization  methods  are  presented.  First, 
the  Korte  nonlinear  least-squares  (LS)  method  is  given.  This  is 
followed  by  two  techniques  of  response  surface  methodology  (RSM) , 
those  being  the  steepest  descent  method  and  the  second-order 
search. 

2.8.1  The  Korte  Non-Linear  Least-Squares  (LS)  Method.  The 

optimization  scheme  used  by  Korte  [8,9]  was  based  upon  the  research 
of  Huddleston  [29] .  The  Korte  solver  optimizes  aerodynamic  design 
through  the  use  of  a  nonlinear  least-squares  optimization 
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formulation.  Basically,  the  nonlinear  optimization  problem 
attempts  to  minimize  an  objective  function  through  determination  of 
design  parameters. 

In  this  method,  the  objective  function,  (P) ,  is  dependent  on 
a  set  of  design  parameters,  (a) ,  and  is  constructed  of  a  number  of 
functions,  (p^)  ,  in  the  nonlinear  least  square  form 

P(*)  =  £p/(*).  (90) 

1=1 

In  vector  form, 

P(a)  -  p  '(a)  p(a)  (91) 

where 

P  =  (Pi.  Pa.  . .  Pj^  (92) 

and 

a  =  (a^,  a^,  . .  a^)^.  (93) 

Thus,  there  are  n  design  parameters  and  m  total  functions. 

The  components  of  p  are  as  follows:  for  i  =  1  to  mi,  the 
weighted  residual  (or  Pi)  is  the  error  in  axial  Mach  number  at  the 
nozzle  exit  and  is  defined  as 

)  k-j;j-l,m^  (94) 

where  co^  =  4 . 0/Mdesign,*  from  mi+1  to  ma,  the  weighted  residual  is  the 
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error  in  flow  angle,  0,  at  the  nozzle  exit  and  is  defined  as 


-0) 

<(>  '  ^imax,j  ' 


k  -  j*m^;  j  -  Iria^-m^ 


(95) 


where  co,^  =  2.5;  and,  finally,  from  m2+l  to  m,  the  weighted  residual 
is  the  error  in  the  axial  Mach  number  along  the  centerline  and  is 
defined  as 


P.  =  --"...J 


*i,l 


axis  ‘ 


k  -  i*m^}  i  = 


(96) 


where  co_  =  1.0/Af  ,  .  (The  centerline  Mach  number,  Af_.  ,  distribution 
is  obtained  using  equations  derived  and  readily  available  in  [37].) 

The  design  parameters,  a,  must  now  be  defined.  First,  the 
nozzle  contour  is  specified  using  the  following  n  cubic  polynomial 
equations : 


r(x) 


•^2 


0  12 


^4n-4*^4n-3^ 


'♦a3X^ 

II 

>4 

0  ^  X  ^  Xj 

X  -  x-x^; 

Xj  ^  X  ^  Xj 

X„T  ^  X  i  x„ 
n-l  n 

(97) 


These  n  cubic  equations  produce  4n  coefficients  given  that  the  x^ 
are  known.  The  x^  are  the  x-locations  of  the  ciibic  spline  knots 
mentioned  earlier.  By  requiring  continuity  of  the  surface  and  the 
slope  and  the  curvature  at  the  endpoints  of  the  segments,  3(n-l)  of 
the  coefficients  are  specified.  The  remaining  coefficients  are 
specified  by  the  predetermined  inlet  and  exit  radii  and  slopes  (4 
coefficients)  and  by  assigning  wall  slopes  at  the  n-1  interior 
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points.  The  set  of  equations  is  now  a  linear  system  of  equations 
which  can  be  solved  to  determine  all  of  the  coefficient  values. 
The  design  parameters  are  the  wall  slopes  at  each  of  the  interior 
points  (ai,  as,  ag,  ...  ,  .  In  other  words,  these  are  the 
components  which  make  up  the  a  vector. 

Using  the  Newton  method  for  nonlinear  least  squares  (assuming 
this  is  a  small  residual  problem) ,  the  equation  which  must  be 
solved  is 

(98) 


where  k  is  the  iteration  number  for  the  updated  wall  solution,  and 
J  is  the  Jacobian  matrix,  defined  as 


da^ 


J  - 


9Pi 


(99) 


Once  Aa^  is  found,  the  updated  wall  slope  coefficients  are  obtained 
using 


►Aa, 


(100) 


With  the  updated  wall  coefficients,  a  new  wall  shape  is  obtained 
using  Eq.  97.  With  the  new  wall  shape,  another  CFD  solution  must 
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be  obtained.  This  procedure  is  repeated  until  the  wall  shape 
converges;  that  is,  the  change  between  the  kth  iteration  and  the 
k+lth  iteration  falls  within  some  tolerance  level. 

As  noted  by  Korte  et  al.,  the  computer-time  expensive  portion 
of  this  procedure  is  in  obtaining  the  Jacobian  matrix.  Each 
interior  point  must  have  its  slope  changed  in  order  to  generate  the 
columns  of  the  J  matrix.  Each  slope  change  is  done  independently 
of  the  other  interior  points'  slope  changes,  and  for  each  slope 
change  made,  another  CFD  flow  field  solution  must  be  obtained  in 
order  to  generate  a  column  of  the  J  matrix.  To  illustrate,  if 
there  are  100  interior  points,  101  CFD  solutions  are  needed  for  one 
iteration  (100  +  1  initial  solution) .  So  it  would  appear  to  pay  to 
minimize  the  design  parameters  in  an  effort  to  cut  down  on  CFD 
solutions  needed. 

2.8.2  Response  Surface  Methodologies  (RSM)  Techniques.  RSM 

techniques  represent  a  somewhat  different  approach  to  the  solving 
of  an  optimization  problem.  Though  minimization  of  an  objective 
(or  response)  function  through  variation  of  design  parameters  is 
still  performed,  the  method  in  which  the  path  to  the  minimum  is 
taken  is  determined  by  fitting  a  "mathematical  French  curve"  to  a 
smooth  response  surface  [45].  "By  careful  design  and  analysis  of 
experiments,  it  seeks  to  relate  a  response,  or  output  variable  to 
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the  levels  of  a  number  of  predictors,  or  input  variables,  that 
affect  it"  [45] . 

"Response  Surface  Methodology  (RSM)  comprises  a  set  of 
statistical  and  mathematical  techniques  for  empirical  model 
building  and  exploitation  ..."  [72].  It  includes  the  development 
of  a  group  of  experimental  runs  that  will  produce  "adequate  and 
reliable  measurements  of  the  response (s)  of  interest  in  a  region  of 
interest";  the  analysis  of  the  experimental  results  to  detemine  an 
empirical/mathematical  model  which  provides  a  best  fit  to  the 
response  surface  data;  and  the  determination  of  the  values  of  the 
design  parameters  which  provide  the  desired  response  (possibly  a 
maximiam  or  a  minimum)  [72]  .  Some  explanations  and  definitions 
follow. 

In  RSM,  a  response  is  defined  as  "the  output  of  a  system  or 
process  that  occurs  as  a  result  of  (or  in  response  to)  a  set  of 
inputs"  [72].  The  inputs,  (i  =  l,k),  are  put  through  the  system 
or  process  and  result  in  the  output,  z  (Fig.  8) .  In  mathematical 

language,  z-ft^x^/X^, _ ,Xj^)-fiX),  and  even  more  specifically, 

z-f{X,Q),  where  9  ( 0^,  . . . ,  0^)  represents  the  set  of  physical 
parameters  that  help  to  define  the  relationship  of  the  response 
function.  The  latter,  more  specific  form  of  z  is  referred  to  as 
the  "mechanistic  model".  Unfortunately,  the  true  response 
function,  f,  is  generally  not  known  and  may  be  quite  complex  and 
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therefore  must  be  approximated.  Additionally,  0  itself  is 
generally  unknown  and  must  be  approximated  as  well. 

The  response  surface  is  defined  as  "a  geometric  representation 
of  a  response  function"  [72]  .  An  example  of  a  response  surface  is 
presented  in  Fig.  9  for  a  case  in  which  the  response  surface 
parameter,  yield,  is  plotted  as  a  function  of  the  two  design 
parameters  of  reaction  temperature  and  reaction  time  [45] .  The 
associated  contour  plot  is  given  in  Fig.  10. 

One  tests  the  system  with  different  sets  of  value  of  the 
design  parameters  to  use  the  generated  responses  to  aid  in 
developing  an  "approximate  interpolating  function"  to  the  response 
surface  [72] .  This  approximating  function  would  be  given  by 
z  »  g(x^,  ^2,  . . . ,  3^,  Pj,  . . . ,  3^)  =g(x,  and  is  referred  to  as  the 
"empirical  model."  "The  operation  of  the  system  with  the  k  design 
inputs  adjusted  to  some  definite  set  of  levels  (values)  is  referred 
to  as  an  experimental  run"  [72] . 

The  techniques  used  in  RSM  include  designed  experiments,  in 
which  sufficient  and  purposeful  change  is  brought  about  in  the 
design  parameters  to  allow  observation  and  modelling  of  the  changes 
in  the  response  (s) ;  regression  analysis,  in  which  "statistical 
techniques  (are)  used  to  model  the  response  as  a  linear 
combination"  of  the  input/design  variables,  X,  and  their 
interactions,  0;  and  steepest  ascent /descent,  in  which  a  "gradient 
search  technique"  is  used  to  determine  the  location  of  either  the 
next  point  at  which  to  begin  the  next  search  for  the  desired 
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response,  or  the  point  at  which  a  desired  response  is  obtained 
[72]  . 

A  graduating  function  is  defined  as  "an  empirical  model  that 
approximates  the  mechanistic  model  over  a  specified  region  of 
interest"  [72]  .  A  graduating  function  is  used  "to  approximate  the 
true  response  function"  [72] .  This  requires  a  means  of  fitting  the 
empirical  functions  to  actual  experimental  data. 

Assuming  an  adequate  model  has  been  postulated  to  represent 
the  data,  one  must  then  obtain  the  best  estimates  of  the  model 
parameters  (0)  in  the  attempt  to  fit  the  empirical  functions  to  the 
data.  In  RSM,  the  method  of  least  squares  is  used  to  obtain  the 
best  estimate  of  6.  It  should  be  noted  that  z  is  restricted  to 
being  linear  in  the  parameters,  meaning  that 
z  =  0) -0iyi+02^2*“**®pJp'  where  y2/  t  •  •  -  rVp  "denote  known  constants 

that  are  functions  of  the  input  variables  . . . ,  ..."  [72]. 

Additionally,  while  z  is  restricted  to  being  linear  in  the  y,  0,  ^ 
parameters,  it  is  not  restricted  to  being  linear  in  the  X  input 
variables . 

In  the  designed  experiments  phase,  one  performs  a  "set  of  n 
experiments  on  the  process,  observing  the  response  at  each  of  n 

sets  of  experimental  conditions  denoted  by  _ ,X^,  obtaining 

the  vector  of  responses  z- (z^,  z^, _ [”^2].  In  the  subsequent 

regression  analysis  phase,  the  error  between  the  observed  response 
and  that  calculated  by  the  mathematical  model  is  represented  by 
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Assuming  again  that  the  response  function  is  linear  in  the 
parameters  mentioned  above,  the  mathematical  model  can  be 
represented  by 

z-Ye*e  (102) 
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Least  squares  is  used  to  find  that  value  for  0  which  "minimizes  the 
sum  of  squares  function  defined  as  the  sum  of  the  squared  model 
error  and  given  by"  [72] 

5(0)  =  (z-re)^(z-re)  .  (104) 

The  desired  value  for  0  is  then  obtained  by  setting  ds  (Q)  /dQ-0 
which  produces 

(r'r)0  =  r'z  (105) 

which  represents  a  set  of  p  equations  which  are  normal  and  which 
can  be  solved  for  0  given  that  r'r  is  non-singular.  Having  solved 
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for  0  it  is  now  possible  to  use  steepest  as cent /descent  methods  (as 
described  below)  to  further  the  search  for  the  desired  response 
from  the  system. 

It  is  important  at  this  juncture  to  take  note  of  a  point  of 
clarification  for  the  remainder  of  this  paper.  Both  the  first- 
order  and  second-order  response  surface  methods  described  below  are 
steepest  descent  methods.  However,  the  first-order  method  will  be 
referred  to  as  the  steepest  descent  method,  and  the  second-order 
method  will  be  referred  to  as  the  second-order  search  method. 

In  the  steepest  descent  method  and  the  second-order  search 
method,  though  not  absolutely  necessary,  a  coordinate 
transformation  may  be  performed  in  order  to  facilitate  the  search 
procedure.  As  further  illustration,  assume  for  the  moment  that 
only  two  generic  design  parameters  are  used  and  are  represented  by 
and  .  This  transformation  is  made  to  take  place  by  first 
selecting  a  point,  designated  by  (x, ,  x, )  =  (X,  ,  x,  )  which  is  to 

be  the  center  of  the  search.  Next,  a  delta  to  x^  and  X^, 
designated  by  dx^  and  dX2,  is  chosen.  Then  the  coordinates  for  the 
new  system  are  obtained  from 


dX, 


(106) 


and 
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(107) 


These  new  variables  are  referred  to  as  the  coded  input  variables. 
The  centerpoint  then  becomes  (x^,  x^)  <=>  (0,  0).  For  the  steepest 
descent  method,  the  additional  points  consist  of  the  following 
pairs  of  {X^,X2)  :  (1,1),  (1,-1),  (-1,-1),  and  (-1,1)  (Fig.  11). 

For  the  second-order  search  method,  in  addition  to  the  points  used 
in  the  steepest  descent  method,  the  additional  points  consist  of 
(0,/2),  (0,-v^),  (^/2,0),  (-/2,0)  (Fig.  12). 

With  the  coded  variables  in  use,  Eq.  102  is  now  represented  by 

*  =  (108) 

where  the  form  of  x  is  given  in  the  following  two  sections.  Once 
again  when  minimizing  the  sum  of  the  squares  function,  the 
resulting  equation, 

{X^X)Q,-X*z  (109) 

can  be  solved  for  ^  given  that  x’x  is  non-singular. 

Though  optimization  methods  used  by  other  researchers  have 
concentrated  on  several  design  parameters,  it  was  hoped  that  a 
reasonably  accurate  engineering  solution  to  the  optimized 
hypersonic  nozzle  could  be  obtained  with  the  use  of  only  two  design 
parameters.  By  using  only  two  design  parameters,  it  would 
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hopefully  be  possible  to  use  RSM  in  a  straightforward  manner  to 
obtain  a  minimum  to  the  objective  function,  and  hence,  find  an 
optimized  nozzle  shape.  The  two  design  parameters  chosen  for  this 
research  consisted  of  the  attachment  angle,  ©attach/  arid  the  exit 
angle,  ©g^it  (Fig.  13)  .  It  was  felt  that  the  necessary  nozzle  wall 
shape  to  produce  uniform  nozzle  exit  flow  at  the  design  conditions 
could  be  obtained  with  just  these  two  parameters,  given  that  the 
nozzle  length  and  the  radius  of  the  attachment  circle,  r^irc/  were 
prescribed  by  the  designer  (Fig.  13) . 

2. 8. 2.1  The  Steepest  Descent  Method.  The  steepest  descent 
method  is  a  first-order  method  which  determines  the  direction  of 
largest  change  based  on  the  sampling  of  the  points  mentioned  above. 
Once  the  response  is  determined  for  each  of  those  points,  a  first- 
order  mathematical  approximation  to  the  local  response  surface  is 
obtained,  then  the  direction  of  largest  gradient  is  determined  and 
a  search  is  begun  in  that  direction  to  find  a  new  minimum. 

In  order  to  characterize  the  local  surface  with  a  linear 
equation,  a  first-order  least  squares  bivariate  regression  equation 
of  the  form 

(110) 

must  be  obtained.  In  order  to  do  so,  the  coefficients  bo,  bi,  and 
b2  must  be  determined.  Recall  that  merely  an  artifice 
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used  to  approximate  the  true  local  surface  at  any  point, 
represented  by  F.  In  order  to  best  approximate  the  surface  and 
solve  for  bo,  bi,  and  b2,  the  equation 

X*fb-X^z  (111) 


must  be  solved  for  b,  where 

b-  (bg,  b^, 

‘1 
'  1 

X  -  ; 

.1 


Z  -  (z^f  Z^f  . . . ,  z^)  ^ 


(112) 


(113) 


(114) 


and  N  is  the  total  number  of  points  being  sampled  and  is  the 
response  value  associated  with  point  (x,  ,  x.  )  [45]  .  Note  that  b 
is  synonymous  with  ^  in  Eq.  109  above.  Another  form  of  Eq.  Ill  is 
given  by 

Ab-Z,  (115) 


where 


A  -  X*X 


(116) 
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and 


Z  -  X^z. 

The  matrix  A, a  3x3  matrix  represented  by 


f 

\ 

^11 

^12 

^13 

A  = 

^21 

^22 

^23 

k  ^31 

^32 

^33; 

has  the  following  components: 


For  the  steepest  descent  method  described  above 
the  values  for  N  and  the  prespecified  values 
matrix  a  simply  becomes 


A  = 


0  4  0 
.0  0  4, 


Once  the  matrix  a  is  obtained,  in  order  to 
necessary  to  obtain  A'^  and  also  the  vector  z  (= 
is  defined  by 


yN 


E-l 


(117) 


(118) 


(119) 

I  is  5.  Inputting 
for  (X.  ,  X,  ) ,  the 

(120) 

solve  for  b  it  is 
[z^,  z^,  z^]"^)  which 

(121) 
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In  this  particular  case,  since  the  x,  and  x,  don't  change,  the  A'^ 
matrix  remains  constant  and  is  represented  by 


0.2  0  0 
0  0.25  0 

0  0  0.25 


(122) 


At  this  point,  it  is  possible  to  calculate  the  coefficients  in  Eq. 
110  above,  which  are  obtained  from 

h-A-^Z.  (123) 


Once  h  has  been  computed,  the  direction  of  steepest  descent 
has  been  found  because  the  gradient  vector, 

h' -  (124) 


has  been  determined  [45] .  In  searching  for  a  new  minimum  if  one 
wishes  to  move  a  specific  distance  g  in  the  steepest  descent 
direction,  one  need  only  specify  g,  from  which  the  x^  and  can  be 
calculated  from  the  equations 


sign  (ijj) 


(125) 


and 


X, 


(126) 
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where 


sign[b^) 


Jb^^O 
Jb^>0  * 


(127) 


The  determination  of  and  is  obtained  from  the  fact  that  the 
magnitude  of  the  coefficient  vector  at  any  point  g  along  the 
direction  of  steepest  descent  is  defined  by 


(128) 


and  the  ratio  of  X^/x^  is  the  same  as  the  ratio  of  since  the 
point  lies  along  the  direction  of  steepest  descent  [45] .  Note  that 
sign(bi)  (Eq.  125)  is  important  since  the  sign  of  b^  and  b2  may  be 
lost  in  the  ratio  bs/bi.  Also,  note  that  if  bi  is  negative,  then 
the  direction  of  steepest  descent  requires  that  x^  be  positive  in 
order  to  move  in  the  correct  direction.  The  same  holds  true  for  b2 
and  x^ . 

Now  that  the  direction  of  steepest  descent  has  been  found,  one 
must  sample  points  in  that  direction  until  an  inflection  point  has 
been  detected.  By  sampling  at  least  3  points  in  the  direction  of 
steepest  descent,  it  should  be  possible  to  determine  if  there  is 
indeed  a  new  minimum  in  this  direction  (Fig.  14)  .  If  there  is  a 
new  minimum  in  this  direction  and  the  sampling  is  fine  enough,  then 
a  good  second-order  curve  fit  to  the  response  surface  in  the 
direction  of  steepest  descent  should  be  obtainable  in  trying  to 
find  the  new  minimum. 
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For  example,  if  three  points  were  used  (gi,  gj,  and  g^,  where 
gi  <  g2  <  ga)  with  corresponding  response  functions  (Fig.  14)  )  , 
and,  assuming  a  new  minimum  occurred  at  point  gj,  then  it  should  be 
possible  to  obtain  a  least-squares  second-order  fit  of  the  along 
the  direction  of  steepest  descent  to  the  g^.  In  other  words,  an 
equation  of  the  form 

bl^  *  bl^g  *  bl^g^  (129) 

should  be  obtainable,  where  the  bl  vector,  (bio,  bl^,  bl2)'^,  is 
represented  by 

bl-Al-^zi.  (130) 


The  development  of  these  equations  can  easily  be  seen  by  replacing 
with  g  and  with  g^  in  Eq.  110.  The  components  of  the  A1 
matrix  are  then 


11 

=  N1 

^^12 

3^13 

21 

ali2 

^^22 

=  ali3 

^^23 

31 

=  ^^13 

^^32 

=  ^^23 

3^33 

=  (?.)• 


and  the  zi  (=  [zi^,  zi^,  vector  is  given  by 


Zl^  -  52^1  Zi  21, 


-  51^1  9-z 


Zl. 


- 


?;zi. 


(131) 


(132) 


where  N1  is  the  number  of  points  being  used  to  attempt  a  fit;  in 
this  case,  N1  =  3. 

Following  the  determination  of  the  bl  vector,  since  fl(g)fit  is 
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a  quadratic  equation,  by  taking  the  derivative  of  Eq.  129  with 
respect  to  g  and  setting  that  derivative  equal  to  0  to  find  the 
minimum,  the  following  equation  is  obtained: 

=  <133) 

which,  when  substituted  into  Eq.  129,  produces  the  curve  fit 
minimiam  of  fl(g)fit.  In  reality,  whether  or  not  a  new  minimum 
exists  at  that  value  of  g  is  very  much  dependent  on  how  well  Eq. 
129  has  characterized  the  response  surface.  The  only  way  to 
determine  the  quality  of  the  response  surface  characterization  is 
to  use  that  value  of  g  to  obtain  new  values  of  and  x^,  and  to 
obtain  the  new  response  function  there,  this  to  compare  with  the 
other  points  along  the  direction  of  steepest  descent. 

Assuming  a  new  minimum  is  found  in  the  direction  of  steepest 
descent,  the  associated  point  becomes  the  center  of  the  next 
stencil  and  presumably  a  smaller  stencil  in  non- transformed  space 
is  used  in  the  next  iteration  of  the  steepest  descent  method  (Fig. 
15)  .  The  rationale  for  the  next  stencil  to  be  smaller  is  to 
provide  a  more  localized  surface  to  be  used  in  providing  the 
requisite  data  for  solving  Eq.  110,  and  to  thus  provide  a  more 
accurate  representation  of  the  local  surface. 

If,  alternatively,  no  new  minimum  is  found,  then  a  number  of 
choices  might  be  made.  First,  a  second-order  search  can  be 
performed  since  5  of  the  9  sampling  points  needed  have  already  had 
their  response  functions  determined  (Fig.  12) .  Second,  a  smaller 
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stencil  can  be  used  in  non- transformed  space  (Fig.  16) .  Third,  the 
minimum  of  the  5  stencil  points  can  be  used  as  the  new  center, 
assuming  the  minimiim  doesn't  occur  at  the  old  center  (Fig.  17)  . 
And  fourth,  some  combination  of  those  three  choices  can  be 
attempted. 

As  long  as  the  surface  can  be  characterized  to  a  good 
approximation  by  a  first-order  fit,  the  method  of  steepest  descent 
works  well.  However,  in  areas  where  strong  second-order  effects 
are  evident  it  is  likely  that  the  first-order  empirical  model  of 
the  surface  will  fail  to  find  a  minimum  and  the  second-order  search 
method  will  need  to  be  utilized  [45] . 


2. 8. 2. 2  The  Second-Order  Search  Method.  The  second-order 
method  begins  much  like  the  steepest  descent  method  in  that  a 
number  of  points  including  a  centerpoint  need  be  sampled  in  order 
to  obtain  some  form  of  characterization  of  the  response  surface. 
The  9  points  mentioned  above  are  sampled  and  each  provides  a 
response  value  (Fig.  18) .  In  order  to  approximate  the  surface  with 
a  second-order  polynomial,  a  second-order  least  squares  bivariate 
regression  equation  of  the  form 


f2[x^,x^)  . 


■b^x^. 


*^22(^2) 


(134) 


where  the  b  vector  (  =  [Jb^,  Jb^,  b^,  b^^,  b^^,  is  obtained  from 

A2b  -  Z.  (135) 
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A2  is  a  (6x6)  matrix,  and  both  b  and  Z  are  (1x6)  vectors.  The 
development  of  &2  and  z  is  very  similar  to  that  of  A  and  z  in  Sec. 
2. 8. 2.1,  Eqs.  116  and  117  respectively.  The  difference  between  the 
two  methods  starts  in  the  z  matrix,  which  is  now  given  by 


(136) 


where  N  is  the  total  number  of  points  being  sampled.  For  the 
second-order  search  method,  N  is  fixed  at  9. 

For  a  general  distribution  of  points,  A2  is  represented  by 


a2^=lf 

^21^12 

^22^14 

^23=^14 

^31^13 

^32"^  23 

^33=^15 

^34=^26 

^34=^23 

(137) 

^41^14 

^42=^24 

^43^26 

^44=i:«jq; 

^44=i:^i^i,\ 

^52=^25 

^33=^35 

^54=^45 

^54“53i=l 

^61=^16 

^«2-'^26 

^«3=^3« 

^45-^56 

^44=^45 

and  the  z  vector  is  represented  by 
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(138) 


Due  to  a  well-chosen  distribution  of  the  sampling  points  (Fig.  12), 
the  A2  matrix  simply  becomes 


A2  = 


9 

0 

0 

8 

8 

0 


0 

8 

0 

0 

0 

0 


0  8  8 
0  0  0 
8  0  0 
0  12  4 
0  4  12 
0  0  0 


0 

0 

0 

0 

0 

4 


(139) 


and  it  can  readily  be  seen  that  straightforward  solutions  can  be 
determined  for  bi,  ■b2  and  bi2»  These  solutions  are  given  by 


Jb 


1 


(140) 


For  the  remaining  3  equations,  Gaussian  elimination  is  used  to 
provide  the  following  results: 


(141) 
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(142) 


and 


Z^-8b^,-8b22 


(143) 


Once  b  has  been  computed,  the  surface  may  be  characterized  at 
the  stationary  point  (the  minimum)  by  examining  the  second-order 
partial  derivatives  of  Eq.  134.  The  stationary  point,  if  one 
exists,  is  the  point  at  which  d{f2)/dX^-0  for  all  X^.  By  taking 
the  second-partial  derivatives  of  f2,  the  Hessian  matrix.  He, 
represented  by 


He 


^12  ^^22 


(144) 


is  obtained.  Then,  taking  its  inverse,  He"^,  and  multiplying  by  the 
V  vector,  represented  by  (-bi,-b2)'^,  the  new  stationary  point  is 
obtained  as  follows  [45] 

f  -  \ 


\  2,/ 


He'^V. 


(145) 


At  this  point  in  the  calculations,  a  new  theoretical  minimum  point 
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has  been  established  at  (x  ,  x  )  ,  and  the  only  remaining  step  in 

^3  ^3 

this  iteration  is  to  obtain  the  response  value  at  that  point.  If 
the  stationary  point  provides  a  new  minimum,  then  the  search  is 
begun  anew  with  the  stationary  point  at  the  center  of  the  sampling 
space  and  presiamably  a  smaller  stencil  in  non- transformed 
coordinates  to  try  to  zero  in  on  the  global  minimum  (Fig.  19) . 

Should  the  stationary  point  not  provide  a  new  minimum,  then  a 
number  of  approaches  exist.  First,  one  might  try  a  smaller  stencil 
in  non- transformed  space  (Fig.  20) .  Second,  one  might  choose  the 
response  function  minimum  of  the  nine  sampled  points  in  the  last 
iteration,  assuming  it  doesn't  occur  at  the  old  center,  as  the  new 
center  and  begin  a  new  search  (Fig.  21) .  Third,  one  might  revert 
to  a  first-order  search  and  use  the  steepest  descent  method  to 
perform  the  next  phase  of  the  search.  And  again,  one  might  use 
some  combination  of  the  three  choices  already  presented. 
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3.  Computational  Techniques 


This  chapter  details  the  methods  used  to  convert  the  governing 
equations  into  discretized  form  for  CFD  solution.  The  discrete 
form  of  the  NS  equations  in  perfect  gas  form  is  presented,  followed 
by  a  short  discussion  on  the  discrete  form  of  the  nonequilibrium 
chemical  reaction  (NECR)  NS  equations.  Subsequently,  dis¬ 
cretization  of  the  perfect  gas  PNS  equations,  the  Baldwin-Lomax  (B- 
L)  turbulence  model  and  the  Yang-Shih  (Y-S)  k-e  model  are  given. 
The  boundary  conditions  are  presented  at  the  end  of  each  of  the 
sections  mentioned. 


3. 1  Discretization  of  the  Perfect  Gas  NS  Equations  and  Boundary 
Conditions 

Temporarily  ignoring  the  viscous  and  source  terms,  the 
governing  NS  equation  (Eq.  2)  takes  on  the  form  in  transformed 
coordinates 


Act"  , 


—  ( F  "  f""  " 

^UV2,J  ^1-1/2, J 


} 

An 


'  n 


-G 


-  n 

^i,j-l/2 


) 


(146) 


where 
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(147) 


*^\.U2j\.njKu2,,^  (C;2,r^W2,,)  ] 

The  last  term  on  the  right  hand  side  of  Eq.  147  is  the  dissipative 
term  which  obviates  the  necessity  of  additional  numerical  damping 
[43]  . 


Once  a  solution  is  obtained  to  Eq.  146,  it  is  necessary  to 
transform  Acj”  .  to  Cartesian  coordinates  before  performing  the  time 

•1#  J 

integration.  The  transformation  takes  the  form 


-- 


An°  ,j 

1,1 


(148) 


The  two-stage  Runge-Kutta  [42]  time  integration  takes  the  form 


(149) 


for  the  predictor  phase  and 


n  n*l  —  -lyj  -lyj 
i,1  ~  o 


(150) 


for  the  corrector  phase. 

The  factors  in  the  last  term  of  Eq.  147  represent  A,  the 
linearized  fo2nn  of  the  computational  flux  Jacobian,  and  consist  of 
the  following:  the  matrix  defined  as 
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The  total  enthalpy,  is  defined  as 

The  matrix  is  defined  as 
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where 


and 


b  (u^*v^) 


(156) 


(157) 


(158) 


The  matrix  A^,  whose  only  non-zero  elements  are  on  the  main 
diagonal  and  are  the  eigenvalues  of  A  is  defined  as 


A 


A 


0  0  0 

0  0  E  u+E  v*k^c  0 

0  0  0  E  u*E  V 


(159) 


(Note  that  c  is  the  speed  of  sound  in  Eqs.  151,  156,  157,  and  159) . 

At  certain  points  where  the  flow  goes  sonic  or  enters  a 
stagnation  region,  the  eigenvalues  may  be  near  zero.  Due  to 
numerical  roundoff,  the  result  may  be  a  violation  of  the  entropy 
condition:  the  scheme  may  produce  the  nonphysical  solution  of  an 
expansion  shock,  a  decrease  in  entropy  [7]  .  As  a  result,  an 
entropy  cutoff  is  used  to  ensure  that  the  eigenvalues  never 
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decrease  below  a  certain  level,  defined  as 
resultantly  defined  as 


The  eigenvalues  are 


X  - 


26. 


|X|3|X| 


|X|<5^ 


(160) 


where  6^,  takes  one  of  two  forms.  The  first  form  is  used  in  the 
body-normal  direction  and  is  represented  by 


5 


X 


-  SJ'’- 


|u  -vpi.lOv^Mvnl) 


(161) 


where  5  is  a  constant  value.  For  this  application,  after 
discussions  with  Gaitonde,  the  value  was  made  to  be  0.01  to 
minimize  the  additional  dissipation  associated  with  this  cutoff 
[57]  .  The  second  form,  the  anisotropic  form  of  5^,  is  used  in  the 
streamwise  direction  to  prevent  excess  dissipation  due  to  large 
grid  cell  aspect  ratios  and  is  given  by  [7] 


1  + 

2/3 

f 

U'^’J 

(162) 


where  =  |u 'Vkl  +  c|Vk| .  Note  that  only  the  two  eigenvalues 
associated  with  the  speed  of  sound  were  cut  off  in  each  direction. 
The  eigenvalue  of  the  contravariant  velocity  was  not  altered,  per 
[57].  Similar  formulations  are  used  to  calculate  s’"  i,  where  ri 
replaces  ^  in  the  eigenvector  (and  its  inverse)  and  eigenvalue 
matrices . 
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Global  second-order  spatial  accuracy  is  achieved  through  the 
MUSCL  (Monotonic  Upstream  Schemes  for  Conservation  Laws) 
formulation  in  combination  with  the  minmod  slope  limiter  [7]  .  The 
minmod  limiter  reverts  the  solution  to  first  order  at  extrema  (such 
as  shocks)  contained  within  the  flow,  thus  preventing  spurious 
oscillations.  As  illustrated  in  Fig.  22,  the  unconserved 
variables  are  extrapolated  from  the  cell  center  to  the  cell 
interface  with  the  MUSCL  formulation  of  van  Leer  [12]  in 
conjunction  with  the  minmod  limiter  [13],  in  order  to  obtain  the 
conserved  variables  immediately  to  the  left  and  right  of  the  cell 
interface.  These  are  then  used  to  calculate  the  flux  vectors 
immediately  to  the  left  and  right  of  an  interface  (giving  the  first 
four  terms  on  the  right  hand  side  of  Eq.  147)  .  Subsequently,  Roe¬ 
averaging  is  used  to  compute  the  last  term  in  Eq.  147.  Due  to 
limitation  of  space,  some  details  have  been  omitted;  for  an 
excellent  description  of  the  Roe-scheme  as  pertains  to  perfect  gas 
flows,  the  reader  is  referred  to  the  work  of  Morrison  [39] . 

As  for  the  remainder  of  the  terms  in  Eq.  13,  the  viscous  terms 
are  handled  in  a  central-differencing  manner  and  the  source  terms 
are  always  lagged  by  the  values  from  the  previous  half  time-step. 


3.1.1  Inlet  Boundary  Conditions.  The  inlet  boundary 
conditions  were  taken  directly  from  the  work  of  Korte  et  al.  [9], 
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as  it  was  felt  this  was  an  excellent  source  for  nozzle  subsonic 


inlet  conditions.  The  first  two  columns  of  cells  at  the  inflow 
(Fig.  23)  must  be  computed  using  these  conditions  in  order  to 
maintain  a  second-order  boundary  condition.  The  assumptions  made 
are  that  the  streamwise  mass  flow  rate  is  conserved  and  the  normal 
mass  flow  rate  remains  constant.  The  other  two  inlet  conditions 
are  made  by  specifying  the  reservoir  stagnation  enthalpy  and 
pressure  (h^  and  p^)  and  accounting  for  compressibility  in  p^.  For 
the  testbed  nozzle,  h^  and  p^  were  specified  as  502,000  J/kg  and 
12,410,563  N/m^,  as  recommended  to  the  author  in  [63].  The  mass 
flow  rate  conditions  are  represented  by 


{pu)i^^ 


=  (pa) 


Area,  . 


(163) 


and 


'  i.i 


while  the  second  two  conditions  are  given  by 


Y-l  P,,^  2 


(164) 


(165) 


and 
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(166) 


(pu)]^ 


]/ 
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where 


P' 


T 


(167) 


and 


P.  -  (1. 

Pt  2 


(168) 


The  second  two  conditions,  Eqs.  165  and  166  combine  to  produce 


h^Q 


2 


ipu)lj*{pv)i:l  ^  ^ 
2(Y-1) 


(169) 


In  solving  the  inlet  boundary  conditions,  Eq.  169  ,  a  quadratic, 
must  be  solved  for  p,  .  using  Eqs.  163,  164,  165,  and  167.  Eqs.  167 
and  168  account  for  compressibility  and  are  generated  using  the 
last  computed  values  (i.e.,  from  the  previous  half  time-step)  of  p, 
u,  V,  and  M.  Once  the  density  is  calculated  using  the  quadratic 
above,  then  the  pressure  can  be  calculated  from  Eq.  166.  Following 
that,  the  variables,  u,  v,  ei,  and  E  are  easily  obtained  with 


u 


i.i 


(170) 
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Pi.l^ 
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and. 
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p,  .  [e.  *  —  (uf  .*vf  . 


)] 
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At  the  end  of  each  half  time-step  after  the  computation  of  the 
interior  of  the  flow  is  complete,  the  inflow  boundary  conditions 
are  enforced,  beginning  with  the  second  column  of  cells  {i=2)  (Fig. 
23) .  After  that  column  is  completely  computed,  the  first  column  of 
cells  (i=l)  is  computed. 


3.1.2  Nozzle  Wall  Boundary  Conditions.  The  boundary 
conditions  for  the  wall  (Fig.  24)  are  as  follows:  no-slip,  a  zero 
first-order  pressure  gradient  in  the  body-normal  direction,  p,  and 
a  specified  wall  temperature  which  had  the  following  distribution: 
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(174) 


^  -  J  500  X:s0. 2032m 

^  “  [509. 51-46. six  x>0. 2032m 

This  wall  temperature  distribution  was  provided  to  the  author  [63]  . 

With  those  three  boundary  conditions  and  the  perfect  gas 
assumption  all  of  the  variables  at  the  wall  are  easily  calculated. 
The  velocity  components  for  the  ghost  cells  (the  grid  cells  in  the 
nozzle  wall,  used  to  enforce  boundary  conditions,  j  =  jl)  (Fig.  24) 
are  obtained  from  the  no-slip  condition  and  are  given  by 


u 


-u 


-  -V 


(175) 


while  the  zero  first-order  pressure  gradient  results  in 


Pi.il  -  Pi.n-i  •  (176) 

With  the  temperature  at  the  wall  specified  by  T„aii,  a  simple  average 
produces  the  internal  energy  of  the  ghost  cell  as 

~  •  (177) 

With  those  four  variables  computed,  the  density  is  easily 
calculated  from 
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Pl,jl 
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(178) 


e . 

^i,n 

and  the  total  energy  is  obtained  using  Eg.  173. 

The  inviscid  flux  condition  at  the  wall  is  due  to  pressure 
only.  The  rationale  behind  this  can  be  clearly  seen  in  Eqs.  4  and 
5.  In  the  inviscid  portion  of  these  two  equations,  no-slip  drives 
the  velocity  components  at  the  wall  to  be  zero,  eliminating  all 
contributions  except  the  pressure. 

3.1.3  Centerline  Boundary  Conditions.  Reflection  is  used  as 
the  boundary  condition  at  the  centerline  (Fig.  25)  .  This  means 
that  all  of  the  flow  variables  for  the  ghost  cells  (those  on  the 
other  side  of  the  centerline,  j=l)  are  exactly  equal  from  one  side 
of  the  centerline  to  the  other  except  for  the  v  velocity  component. 
The  V  component  takes  on  the  negative  value  of  that  of  its  mirrored 
cell.  Only  one  row  of  cells  is  used  on  the  other  side  of  the 
centerline.  The  flux  condition  at  the  centerline  is  assumed  to  be 
zero  since  the  formulation  is  finite  volume  and  presumably  no  flux 
can  occur  through  the  centerline. 

3.1.4  Nozzle  Exit  Boundary  Conditions.  The  boundary 
conditions  at  the  exit  (Fig.  26)  are  to  assume  a  second-order  zero 
gradient  in  the  ^  direction  for  the  conserved  variables,  as 
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recommended  to  the  author  by  Korte  in  one  of  many  discussions  [53] . 
This  assumption  results  in 


^11,  j  ^ 


(179) 


3.2  Discretization  of  the  HIECR  NS  Equations  and  Boundary  Conditions 

Similar  to  the  perfect  gas  case,  formulations  for  the 
nonequilibrium  gas  case  for  the  Roe-averaged  variables  and  the 
flux-difference  splitting  have  been  developed  and  are  available. 
For  excellent  examples,  the  reader  is  referred  to  [35,  36] . 

3.2.1  Boundary  Conditions.  The  boundary  conditions  for  the 
inflow  are  handled  exactly  as  they  are  for  the  perfect  gas  case  in 
Section  3.1.1,  but  with  the  assumption  that  the  air,  though  at  high 
temperature  and  pressure,  is  at  equilibrium  and  of  a  fixed 
composition.  Thus,  the  vibrational  temperature  is  assumed  to  be 
equal  to  the  temperature  of  the  air,  and  the  mass  fractions  of  N2, 
O2,  N,  0,  and  NO  are  assumed  to  be  constant  at  0.767,  0.233,  l.E- 
06,  l.E-06,  and  l.E-06  respectively.  All  other  boundary  conditions 
are  done  exactly  as  in  the  perfect  gas  case  in  Sections  3.1.2- 
3.1.4. 
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3.3  The  Discrete  Form  of  the  PNS  Equations  in  Perfect  Gas  Form 

The  goal  was  to  produce  a  PNS  solver  which  would  match  very 
well  with  the  NS  solver  and  which  would  have  incorporated  within  it 
the  same  turbulence  models  as  the  NS  solver.  The  reason  that  it 
was  important  that  the  two  solvers  match  well  was  that  the  NS 
solver  was  to  provide  the  input  to  the  PNS  solver  (Fig.  27) .  Thus, 
the  same  methodology  was  used  to  minimize  the  differences  when 
transferring  from  NS  solver  to  PNS  solver.  In  other  words,  the  two 
codes  would  have  in  common  the  same  grid  generation  technique,  the 
same  flux-difference-split  calculations  in  the  p-direction,  and  the 
same  turbulence  model  formulations.  The  PNS  solver  is  overall 
second-order  accurate  and  uses  the  same  two-stage  Runge-Kutta 
integration  scheme  as  the  NS  solver  [42] . 

Ignoring  the  viscous,  source,  and  GCL  terms,  the  parts  of  the 
main  governing  equation  dealing  with  convection  and  pressure,  Eq. 
45,  becomes 


An 


(G~"  -G/  ) 


(180) 


where  retains  the  same  definition  as  in  the  NS  formulation, 

^l,j±l/2 

similar  to  Eq.  147.  The  viscous  terms  are  again  centrally- 
differenced.  The  source  terms  are  once  again  lagged  from  the 
previous  plane  (Fig.  28)  .  It  is  most  important  to  retain  the 
differencing  of  the  metrics  in  the  GCL  terms  in  the  same  form  in 
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which  the  flux  differencing  is  computed,  otherwise  instability  to 
the  solution  may  arise  [9] . 

Once  the  F*  flux  terms  have  been  calculated,  the  conserved 
variables  must  be  decoded  from  the  F*  terms  [49]  .  With 

F‘  -  [pu,  puu*(i>p,  puv,  (E*p)u]''^-  [F*,  F^,  F3,  F^']^  (181) 

the  conserved  variables  are  calculated  from 


where 


V 

Without  the  stipulation  that  ^  ^ (x)  only,  the  decoding  would  be 

much  more  complicated. 

The  value  for  Vigneron's  coefficient,  w  (Eg.  49),  must  be 
known  before  a  solution  can  be  found.  In  this  application,  (o  is 
lagged  by  defining  it  as  a  function  of  the  previous  space  marching 
step  solution.  Once  calculated,  the  value  remains  unchanged 
through  both  the  predictor  and  corrector  phases  of  the  solution. 
Any  attempt  to  change  <0  in  the  current  space  marching  step  lead  to 
instability  and  could  not  be  remedied.  It  is  strongly  suspected 
that  this  is  due  to  the  inability  of  the  flux-differencing  to 
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account  for  this  change  since  the  eigenvalues  are  calculated  based 
on  the  conserved  variables  and  not  on  the  F*  components,  per  Korte 
[49]  .  In  other  words,  since  the  subroutine  which  performed  the 
flux-difference  splitting  used  the  computational  flux  Jacobian, 
B  -  dG/dU,  and  not  the  formulation  used  by  Korte,  B  -  dG/dF‘ ,  there 
was  no  accounting  for  a  change  of  o)  within  the  numerical 
dissipation  portion  of  Eq.  180. 

As  to  the  space  marching  step  itself,  given  that  the  inviscid 
upwind  algorithm  has  a  Courant-Friedrichs-Lewy  (CFL)  [48]  linear 
stability  limit  of  1  and  a  viscous  limit  of  1/2  (valid  for  2-D  or 
axisymmetric) ,  then  the  allowable  space  step  is  defined  by  [49] 


Ax  ^ 


_ wall _ 

wall  ^  2p 

M^^JY  ^ ^y vail 


(184) 


where  Ay„aii  is  the  change  in  y  from  the  wall  to  the  first  grid  point 
off  of  the  wall  (Fig.  29)  [49]  .  All  of  the  other  variables  in  the 
above  equation  were  taken  at  the  first  cell-center  point  away  from 
the  wall.  This  step  is  applied  with  a  safety  factor  of  0.9  since 
the  geometry  in  this  application  is  a  simple  one  [49] . 

With  an  adaptive  grid  subroutine  provided  to  the  author  by 
Korte  [53],  the  Ay^a^^  (Fig.  29)  is  computed  dependent  on  whether  the 
flow  is  laminar  or  turbulent.  Within  the  siibroutine,  the  value  for 
y  at  the  jth  point  is  given  by 
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where 


^ci  = 


u„-j) 


C  =  ;  3  =  1.02 , 

(j„-i)  "" 


(186) 


The  stretching  parameter  is  allowed  to  vary  in  order  to  meet  the 
^Ywaii  specification. 

Should  the  flow  be  laminar,  the  number  of  points  inside  the 
boundary  layer  is  specified  by  user  input  which  in  turn  aids  in 
computing  the  Ay„aii.  This  value  is  currently  set  to  20  points 
inside  of  the  laminar  boundary  layer.  Essentially,  how  this  works 
is  the  subroutine  in  this  case  takes  l/20th  of  the  Mach  number  at 
the  centerline  of  the  nozzle  and  drives  6  such  that  the  first 
point  away  from  the  wall  matches  that  value  of  Mach  number  (Fig. 
30) .  Though  this  uses  a  somewhat  simple  linear  approximation  to 
the  Mach  number  distribution  in  the  boundary  layer,  it  appears  to 
provide  excellent  results  in  obtaining  the  specified  number  of 
points  within  the  boundary  layer  [53] . 

Alternately,  should  the  flow  be  turbulent,  the  value  of 
the  y’  value  (Eq.  54)  of  the  first  point  off  of  the  wall  (Fig.  31), 
is  a  user  input;  the  specified  value  of  in  turn  aids  in 

computing  for  the  first  point  off  of  the  wall.  In  this  case. 
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the  subroutine  drives  such  that  the  first  point  away  from  the 
wall  matches  the  user-specified  value  of  •  The  required  value 
of  varies  depending  on  how  accurate  a  solution  is  deemed 
necessary  by  the  user. 

3.3.1  Boundary  Conditions.  The  boundary  conditions  for  the 
PNS  solver  were  developed  or  taken  from  previously  mentioned 
boundary  conditions.  The  NS  solver  produces  two  adjacent  columns 
of  cells  which  are  the  initial  conditions  for  the  PNS  code  (Fig. 
27)  .  The  most  critical  point  here  is  for  the  NS  solver  to  give  an 
initial  starting  plane  which  does  not  violate  as  will  be 
further  discussed  [9]  .  This  is  done  by  specifying  for  the  NS 
solver  the  value  of  the  safety  factor,  a  (Eg.  49)  ,  letting  the 
solver  backmarch  from  the  exit  through  the  flow  field  solution  in 
search  of  the  plane  at  which  the  violation  occurs,  and  taking  the 
next  column  downstream  as  the  initial  plane  for  the  PNS  solver. 
The  wall  and  centerline  conditions  are  handled  exactly  as  they  are 
in  Sections  3.1.2  and  3.1.3. 

No  outflow  conditions  are  required  for  a  PNS  solver  in 
supersonic  flow.  The  last  solution  plane  at  the  nozzle  exit 
provides  the  exit  boundary  conditions.  In  essence,  this'  is  a 
zeroth-order  extrapolation  or  boundary  condition. 
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3.4  Discretization  of  the  B-L  Turbulence  Model  and  Boundary 
Conditions 

In  solving  for  the  eddy  viscosity  with  the  B-L  turbulence 
model,  it  is  first  necessary  to  find  the  maximum  vorticity, 

This  almost  always  occurs  at  the  wall,  but  to  ensure  it  is  located, 
a  search  marching  away  from  the  wall  for  a  given  axial  location  is 
undertaken.  With  known,  the  value  for  F(y),  Eq.  58,  at  each 
point  can  be  calculated.  Once  that  is  done,  the  point  at  which 
F(y)  has  a  maximiim  value  determines  F„ax  and  y„ax  for  that  column  of 
cells  (Fig.  32) .  Next,  p  and  p  are  calculated  marching 

^iiUisr  ^oatar 

away  from  the  wall.  Until  the  point  at  which  p  is  greater  than 

^iaiiar 

p  ,  the  eddy  viscosity  is  equal  to  p  ,  otherwise  it  is  equal 

^oatar  ^iaaar 

to  p  (Fig.  33)  . 

^outmr 

In  general,  the  niomber  of  times  the  B-L  turbulence  model  is 
updated  in  the  time-integration  scheme  has  an  influence  on  the 
numerical  solution..  However,  the  final  converged  numerical 
solution  at  steady  state  is  independent  of  the  nuonber  of  times  the 
B-L  model  is  updated  within  the  solver.  The  B-L  turbulence  model 
is  therefore  called  once  every  5  total  time  integrations  (that  is, 
once  every  5  predictor-corrector  cycles)  in  the  NS  codes,  per 
recommendation  by  Shang  [64] .  This  is  done  to  aid  in  preventing 
the  NS  solver  from  reacting  too  quickly  to  immediate  changes  in  the 
eddy  viscosity.  In  the  PNS  code,  the  B-L  turbulence  model  is 
called  after  each  spatial  integration  step  (either  predictor  or 
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corrector) ,  since  there  is  no  temporal  dependence  in  the  PNS 
solution,  and  the  geometry  and  physics  change  the  values  of  C0n,ax/ 
F„ax/  and  y„ax  at  each  spatial  location. 

3.4.1  Boundary  Conditions.  There  are  only  two  boundary 
conditions  which  need  to  be  accounted  for  with  the  B-L  model.  At 
the  wall,  eddy  viscosity  equals  0,  so  the  equation  for  the  eddy 
viscosity  for  the  wall  ghost  cell  is  handled  the  same  as  the 
velocity  components  in  Section  3.1.2.  The  eddy  viscosity  for  the 
centerline  ghost  cell  is  equal  to  that  of  its  mirrored  cell,  just 
as  in  Section  3.1.3,  where  the  values  of  the  scalar  parameters 
remain  unchanged  across  the  centerline  (Fig.  34) . 

3.5  Discretization  and  Boundary  Conditions  of  the  Yang-Shih  k-e 
Model 

The  discretization  of  a  k-e  model  using  Roe's  flux-difference 
splitting  is  readily  available  in  [39] .  The  discretization 
recommended  by  Morrison  is  exactly  the  one  used.  A  couple  of  key 
points  must  be  mentioned  however. 

First,  in  order  to  prevent  negative  values  of  k  and  e,  minimum 
values  must  be  established  beyond  which  pk  and  pe  are  not  allowed 
to  transgress.  However,  care  must  be  taken  in  assigning  these 
values  since  the  eddy  viscosity  is  proportional  to  the  ratio  of 
k^/e.  The  values  of  pk^n  and  therefore  should  have  some  basis 
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in  reality  or  at  least  a  reasonable  value  compared  to  the  ratio. 
In  this  fomulation,  the  values  for  pk„in  and  are  chosen  to  be 
l.OE-24  times  the  values  of  pk  and  pe  at  the  centerline  in  the 
inflow  of  the  nozzle,  per  recommendations  by  White  of  NASA  Langley 
[65]. 

Secondly,  since  k  is  included  in  the  total  energy  equation, 
care  must  be  taken  to  account  for  this  inclusion  in  the  Roe  scheme 
formulation  [39] .  Of  specific  importance  is  the  fact  that  the 
enthalpy  and  the  speed  of  sound  change  as  a  result  of  this 
inclusion.  The  total  enthalpy,  (Eq.  155) ,  is  now  defined  as 
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which  using  Eq.  26  makes  the  speed  of  sound,  c. 


c  r  [(Y-1)  {u^*v^)-k)] 
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Both  of  these  equations  are  used  in  the  formulation  of  the  Roe 
scheme  in  Eqs.  146  and  147. 

The  method  of  solution  for  the  governing  partial  differential 
equation  with  the  added  k-G  equations,  Eq.  84,  is  to  first  solve 
for  the  NS  equations  throughout  the  domain  of  interest  and  to 
update  the  NS  CFD  solution.  In  other  words,  the  new  values  for  the 
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subsequently  enforced.  During  the  preceding  processes,  the  last 
computed  (that  is,  the  last  half  step)  values  of  k  and  Pt  used. 
Then  using  the  previous  half  step  values  for  the  pertinent  flow 
variables,  the  portion  of  Eg.  84  represented  by  Eqs.  61  and  62  are 
solved  throughout  the  flow  field  in  order  to  obtain  the  most 
current  values  of  k,  e,  and  Pt  •  After  enforcement  of  the  k-e 
boiindary  conditions,  this  half  step  (either  predictor  or  corrector) 
is  concluded  and  the  next  half  step  starts. 


3.5.1  Boundary  Conditions.  The  inflow  boundary  conditions  are 
not  quite  as  straightforward  as  they  have  been  for  other  portions 
of  the  code.  For  the  inflow,  the  eddy  viscosity  was  assumed  to  be 
equal  to  the  molecular  viscosity  at  each  point  of  the  inflow. 
Assuming  a  turbulence  intensity  of  =  0.01  [41,59]  where  the 
turbulence  intensity  is  defined  as 


= 


u 


irJ 


(190) 


the  value  for  k^^j  can  be  readily  computed.  Then,  using  the  value 
for  eddy  viscosity  and  k  at  each  point  along  the  inflow,  e  can  be 
calculated  using  Eq.  67  [41]  .  In  the  Y-S  formulation,  since  e 
can’t  be  solved  for  explicitly  in  the  eddy  viscosity  calculation 
(Eq.  67) ,  it  is  necessary  to  perform  a  Newton  iteration  to  obtain 
G. 
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The  boundary  conditions  at  the  wall  are  as  follows:  k  is  zero 


at  the  wall  and  [52],  represented  by 


=  2v 


d/k 
dy  j 


wall 


with  the  simple  average  given  by 
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results  in 
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The  reflection  boundary  condition  was  again  used  at  the  centerline 
as  in  Section  3.1.3,  and  a  zero  second-order  streamwise  gradient 
was  used  at  the  outflow  to  obtain  pk  and  pe  there,  as  in  Section 

3.1,4. 
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4.  Nozzle  Flov  Code  Development/Validatlon 


This  chapter  details  all  of  the  steps  involved  in  taking  a  2-D 
Roe-scheme  NS  code  to  full  axisymmetric,  modifying  the  code  by 
adding  the  two  turbulence  models,  and  changing  the  resultant  code 
over  to  PNS  for  use  in  the  supersonic  section  of  the  nozzle,  along 
with  the  validation  steps  taken  along  the  way.  Results  of  the 
validations  are  presented  immediately  following  the  presentation  of 
validation  case  details.  Justification  as  to  why  procedures  were 
followed  or  not  followed  is  provided,  as  are  any  pitfalls 
encountered  along  the  way. 

Salient  features  of  all  of  the  NS  and  PNS  codes  are  as 
follows:  finite  volume,  capable  of  being  run  2-D  or  axisymmetric, 
using  explicit  integration,  and  globally  second-order  accurate 
through  the  use  of  flux-difference-splitting  and  the  minmod 
limiter.  In  addition,  the  NS  codes  can  utilize  local  time  stepping 
for  a  steady-state  problem  in  order  to  accelerate  convergence  to  a 
solution.  The  nozzle  flow  codes  are  highly  vectorized,  resulting 
in  their  running  much  faster  on  a  Cray  or  similar  computer. 
However,  the  solvers  can  be  run  on  any  UNIX  machines,  including  Sun 
Microsystems.  In  the  harshest  case  (nonequilibrium  chemically 
reacting  (NECR)  flow,  NS  solver,  for  almost  the  entire  length  of 
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a  ten  foot  nozzle) ,  the  nozzle  flow  codes  require  less  than  8  Mw  of 
computer  memory  on  a  Cray,  or  263  words  per  grid  point. 

4.1  2-D  NS  Perfect  Gas  Code  Modified  to  2-D/Full  Axisymmetric 

In  order  to  have  a  2-D/axisyinmetric  perfect  gas  code  with  the 
appropriate  turbulence  models  and  changeable  source  code,  it  was 
determined  that  the  most  efficient  process  would  be  to  modify  an 
existing  local  2-D  code  which  utilized  the  Roe-scheme.  Though  a  2- 
D/full  axisymmetric  NS  code  was  available  locally,  this  code  was 
designed  to  account  for  vibrational  nonequilibrium  and  chemical 
reaction  effects;  in  the  laminar  case,  for  instance,  this  code 
solved  for  eleven  conserved  variables  at  each  grid  point  (Eqs.  2 
and  27)  vs.  four  conserved  variables  at  each  point  for  a  strictly 
perfect  gas  NS  solver  (Eqs.  2  and  3) .  It  was  therefore  felt  that 
this  code  would  spend  much  wasted  computer  time  in  needless 
computations  when  solving  the  perfect  gas  NS  equations.  However, 
since  the  nonequilibrium  chemical  reaction  code  and  the  original  NS 
perfect  gas  2-D  code  were  developed  by  the  same  person,  Gaitonde 
[61,  66],  and  the  foundation  for  the  nonequilibrium  code  was  the  NS 
perfect  gas  2-D  code,  it  was  felt  that  the  2-D  code  could  be 
efficiently  and  quickly  converted  to  being  alternatively  axi¬ 
symmetric. 

The  NS  perfect  gas  2-D  code  was  a  proven  finite  volume  code 
which  incorporated  the  Roe-scheme  [7] .  The  finite  volume  method 
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ensures  that  the  basic  quantities  of  mass,  momentum  and  energy  are 
conserved  at  the  discrete  level  [42].  Roe's  scheme  was  determined 
to  be  desirable  for  a  niomber  of  reasons  [7]  .  First,  since 
upwinding  is  used,  no  additional  numerical  dissipation'  is  needed. 
Second,  many  other  schemes  are  highly  dissipative  in  viscous  flows, 
giving  large  errors  in  heat  transfer  and  other  undesirable  effects. 
Third,  the  method  is  extremely  stable  and  robust,  both  highly 
desirable  qualities.  Fourth,  the  scheme  does  an  excellent  job  of 
capturing  shocks. 

The  form  of  Roe's  flux-difference  split  method  used  in  the 
computer  code  is  MUSCL  (Monotonic  Upstream  Schemes  for  Conservation 
Laws) .  Essentially,  in  standard  flux-differencing  (non-MUSCL) ,  the 
flux  vectors  (F  and  G  in  Eq.  2)  at  the  cell  centers  are 
extrapolated  to  the  cell  interfaces.  In  the  MUSCL  approach,  the 
conserved,  dependent  vector  (U  in  Eq.  2)  or  the  unconserved, 
dependent  vector  (p,  u,  v,  p)  at  the  cell  center  is  extrapolated 
to  the  cell  interface  in  order  to  obtain  the  flux  vectors  at  the 
interface  [11]  . 

The  perfect  gas  computer  code  as  given  to  the  author  was  2-D 
and  hard-wired  to  compute  the  laminar  flow  about  an  airfoil.  The 
code  had  to  be  converted  to  solve  either  axisymmetric  flows  or  2-D 
flows,  as  mentioned  previously,  and  be  modified  to  compute  nozzle 
flow.  Per  the  recommendation  of  Gaitonde  [57]  and  following  the 
example  set  forth  in  [66],  Eq.  2  was  selected  as  the  governing 
equation. 
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Thus,  the  change  in  the  2-D  code  to  axisymmetric  entailed 
changing  the  metrics  to  account  for  variation  with  the  normal 
variable,  y  Y^y/J/  yrix/J/  yriy/J  for  the  axisymmetric  case 
vs.  ^^/J,  ^y/J,  T]^/J,  riy/J  for  the  2-D  case);  adding  axisymmetric 
terms  to  the  viscous  terms  to  account  for  non-Cartesian  coordinates 
being  used  (the  terms  with  a  factor  of  6  in  Eqs.  7  and  9);  and 
adding  the  source  terms  to  the  right-hand  side  of  Eq.  2.  In 
addition,  since  the  code  lacked  the  anisotropic  entropy  cutoff  in 
the  streamwise  direction  (Eq.  162),  this  formulation  was  added  to 
the  code;  recall  this  cutoff  alleviates  large  amounts  of 
dissipation  which  arise  with  large  grid  cell  aspect  ratios. 

In  addition  to  reformulating  the  code  to  solve  for  axi¬ 
symmetric  flows,  it  was  necessary  to  develop  initial  conditions  for 
the  nozzle.  The  nozzle  flow  was  initially  assumed  to  be  uniform 
sonic  flow.  However,  with  the  boundary  conditions  previously 
presented  in  Sections  5.1.1  through  5.1.4,  the  flow  became  subsonic 
in  both  the  convergent  and  divergent  sections.  As  a  result,  a  1-D 
approximation  to  the  Mach  number  was  used  to  speed  up  convergence 
and  to  ensure  that  the  flow  remained  supersonic  in  the  divergent 
section  of  the  nozzle.  With  the  reservoir  pressure  and  temperature 
specified  and  assuming  the  flow  was  sonic  (M  =  1)  throughout  the 
length  of  the  throat  (Fig.  35) ,  then  the  temperature,  pressure,  and 
velocity  components  throughout  the  throat  could  be  readily  obtained 
from  [54] 
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It  was  necessary  to  perform  a  Newton  iteration  to  match  the  area  of 
the  nozzle  with  the  Mach  niomber  when  dealing  with  the 
convergent/divergent  sections  of  the  nozzle  (Fig.  6),  using  [54] 
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In  Eg.  197,  A  is  the  cross-sectional  area  at  the  axial  location 
where  M  is  to  determined,  and  A*  is  the  cross-sectional  area  of  the 
throat.  Once  M  had  been  obtained  for  a  given  axial  location,  Eqs. 
194,  195,  and  196  were  used  to  obtain  the  pertinent  flow  variables. 

At  this  point,  it  was  desired  to  compare  the  reformulated 
perfect  gas  code  with  the  nonequilibrium  code  being  run  in  perfect 
gas  mode.  However,  the  nonequilibrium  code  was  hardwired  to 
compute  external  flows  when  received  [66] .  The  boundary  conditions 
were  therefore  changed  to  account  for  nozzle  flow  (see  Sect. 
3.2.1).  In  addition,  some  other  minor  modifications  were  made  to 
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make  the  code  match  up  better  in  perfect  gas  mode.  These  included 
the  addition  of  Keyes'  Law  of  Viscosity  (Eq.  20)  for  hypersonic 
nozzle  flow  [63];  the  hard-wiring  of  the  gas  constants  to  match 
those  used  in  the  perfect  gas  code,  specifically,  R,  C^,  and  Cp,-  the 
zeroing  out  of  the  vibrational  energy  terms  with  a  factor,  S^ecr/ 
similar  to  the  2-D/axisyinmetric  switch  6,  which  was  0  in  perfect 
gas  mode  and  1  in  nonequilibrium  mode  (in  Eq.  27,  e^  — >  SnecrSv)  > 
and  the  addition  of  perfect  gas  boundary  conditions  (see  Sects. 
3.1.1  -  3.1.4) . 

4.1.1  Validation  of  Perfect  Gas  NS  Solver,  Laminar  Axi- 
symmetric  Case.  To  validate  the  perfect  gas  NS  solver,  the  two  NS 
codes  (perfect  gas  and  nonequilibrium/perfect  gas)  were  run  in 
perfect  gas  mode  on  a  rather  coarse  53x28  (i,j)  mesh  of  the  WL  Mach 
12  nozzle  (Fig.  36)  .  The  flow  field  solutions  obtained  were 
compared  with  the  nonequilibrium  solution  on  the  same  mesh  and  with 
SAIC  CFD  data  mentioned  earlier  [22]  .  These  comparisons  were  also 
made  to  determine  whether  nonequilibriiim  and  chemical  reaction 
effects  would  need  to  be  accounted  for  in  the  final  optimization 
code.  Direct  comparisons  would  also  help  to  debug  the  codes.  The 
case  run  was  laminar,  since  no  turbulence  models  had  been  added  to 
the  codes  at  this  time.  Convergence  was  achieved  when  1.5 
characteristic  time  units  had  been  achieved,  or  when  the  solution's 
total  residual  dropped  6  orders  of  magnitude.  This  was  the  general 
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criteria  for  convergence,  unless  otherwise  specified  below. 
However,  in  all  cases  the  residual  dropped  at  least  2-3  orders  of 
magnitude  in  achieving  a  converged  solution. 

A  characteristic  time  unit  is  represented  by 

(198) 

where 

=  (199) 

and  t  is  the  total  time  over  which  a  solution  has  been  integrated. 
The  variables  and  vary  with  the  case  and  the  mesh  used.  For 
the  laminar  case  on  this  grid,  was  3.048  m  and  was  the  u 

velocity  component  in  the  nozzle  throat  obtained  using  Eqs.  194  and 
196. 


4. 1.1.1  Results  of  Validation  of  Perfect  Gas  NS  Solver, 
Laminar  Axisymmetric  Case.  The  final  results  from  the  validation 
done  with  the  NS  perfect  gas  solver  and  the  NS  nonequilibrium 
solver  in  perfect  gas  mode,  as  well  as  with  the  NS  nonequilibrium 
solver  in  NECR  mode  were  very  encouraging.  Plots  of  the  centerline 
static  pressure  (Fig.  37)  and  temperature  (Fig.  38)  vs.  the  axial 
direction  showed  no  significant  difference  between  the  solutions  in 
perfect  gas  form  and  that  in  NECR  form.  Examining  Figs.  39  and  40, 
plots  of  the  same  flow  variables  show  only  a  mild  difference  in  the 
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throat  region,  presiomably  due  to  nonequilibrium  effects.  Otherwise 
all  three  solutions  converged  to  the  same  result  as  the  flow 
traversed  the  nozzle.  Thus,  though  nonequilibrium  effects  were 
seen  in  the  throat  region,  the  effect  on  the  nozzle  exit  flow  is 
indiscernible;  the  author  thus  concluded  that  in  the  laminar  case 
nonequilibrium  effects  were  not  significant  for  the  nozzle  in 
question. 

Plots  of  the  Mach  number  (Fig.  41)  and  Pitot  pressure  (Fig. 
42)  at  the  nozzle  exit  similarly  reveal  no  discernible  difference 
between  all  of  the  solutions.  In  addition,  when  comparing  with 
similar  data  computed  by  SAIC  [22]  for  the  laminar  case,  the  nozzle 
exit  results  were  very  comparable.  Again,  with  these  results  the 
author  concluded  that  nonequilibrium  effects  were  not  significant 
enough  in  this  case  to  warrant  inclusion  of  such  effects  in  a 
laminar  solution  of  this  nozzle.  In  addition,  the  author  concluded 
that  the  perfect  gas  NS  solver  and  the  nonequilibrium/perfect  gas 
NS  solver  had  been  validated  in  the  laminar  hypersonic  nozzle  case. 

4.2  Turbulence  Models  Incorporated  Into  Perfect  Gas  and  NECR 
Codes 

In  order  to  determine  whether  the  zero-equation  Baldwin-Lomax 
(B-L)  model  or  the  two-equation  Yang-Shih  (Y-S)  k-e  model  would 
best  model  the  turbulence  in  the  nozzle  flow,  it  was  necessary  to 
incorporate  both  into  the  perfect  gas  code  and  the  nonequilibrium 
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code.  In  defining  the  best  model,  the  objective  was  to  ascertain 
which  model  would  provide  the  turbulent  nozzle  flow  which  best 
compared  with  nozzle  experimental  data.  It  was  possible  that  some 
trade-offs  between  accuracy  and  central  processing  unit  (cpu)  time 
might  need  to  be  made,  especially  in  light  of  the  determination 
that  the  solvers  used  in  the  subroutine  be  as  efficient  as 
possible.  The  B-L  model  would  undoubtedly  be  faster,  since  fewer 
computations  were  necessary,  but  its  accuracy  had  yet  to  be 
determined.  The  Y-S  k-e  model  might  prove  more  accurate  but  its 
cost  in  computer  time  was  as  yet  unknown.  These  were  items  which 
needed  to  be  addressed  before  a  choice  could  be  made  as  to  which 
turbulence  model  to  use  throughout  the  optimization. 

The  B-L  model  added  to  the  codes  was  primarily  adapted  from  a 
previously  developed  modification  to  the  same  2-D  perfect  gas  code 
as  initially  used  here  [67]  .  With  a  minor  modification,  this 
subroutine  is  the  one  used  in  the  code  when  the  B-L  model  is  used; 
the  minor  modification  consisted  of  removing  a  portion  of  the 
subroutine  which  dealt  with  recirculating  flow.  It  was  felt  that 
this  portion  of  the  subroutine  was  unnecessary  since  flow 
separation  and  recirculation  were  not  likely  to  occur  in  the  nozzle 
flow  being  researched. 

The  k-G  subroutine  used  in  the  code  was  developed  by  the 
author  and  was  based  on  the  work  and  turbulence  model  of  Yang  and 
Shih  [52] .  In  order  to  simplify  debugging  of  the  code,  a  Steger- 
Warming  (S-W)  scheme  was  initially  used,  and  the  Launder-Spaulding 
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(L-S)  k-G  model,  used  by  Rizzetta  [41,59],  was  the  k-e  model  used. 
The  S-W  scheme  was  used  initially  because  of  its  simplicity  [55] . 
The  L-S  model  was  initially  used  because  the  author  was  provided  a 
copy  of  Dr  Rizzetta 's  Beam-Warming  2-D  solver,  in  which  the  L-S 
turbulence  model  is  used.  Details  of  this  code  are  available  in 
[41] .  Hence,  for  debugging  the  2-D/axisymmetric  perfect  gas  solver 
and  for  comparisons  with  the  solution  obtained  from  the  Rizzetta 
code,  the  S-W  scheme  with  the  L-S  k-e  model  was  initially  used. 
Once  the  debugging  of  the  S-W  scheme  with  the  L-S  k-e  model  was 
completed,  a  switch  was  made  to  the  Roe  scheme,  but  the  use  of  the 
L-S  model  was  maintained.  Finally,  once  all  the  bugs  had  been 
worked  out  of  that  code,  the  Roe  scheme  was  used  with  the  Y-S  k-e 
turbulence  model  in  the  last  revision  of  the  NS  code  with 
turbulence. 

All  of  the  work  to  this  point  had  been  done  on  the  perfect  gas 
code.  Later,  however,  both  turbulence  models  were  included  in  the 
nonequilibrium  code  so  that  the  comparison  of  the  perfect  gas 
solutions  with  the  nonequilibrium  solution  could  be  made  for  the 
nozzle,  as  discussed  in  Section  4.2.3. 

The  turbulence  model  validation  occurred  in  three  phases.  The 
turbulence  models  and  their  associated  NS  codes  were  first 
validated  in  2-D  against  another  CFD  solution.  The  turbulence 
models  and  codes  were  then  validated  in  2-D  against  experimental 
data.  And  finally,  the  turbulence  models  and  codes  were  validated 
in  axisymmetric  form  against  WL  Mach  12  nozzle  data. 
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4.2.1  Turbulence  Model  Validation  Against  a  CFD  Solution. 

Both  turbulence  models  were  compared  to  previously  validated  CFD 
results  from  Rizzetta  [68] .  The  test  case  was  a  flat  plate  at  zero 
angle  of  attack  (a)  which  the  Rizzetta  code  had  been  previously 
used  to  solve.  Adiabatic  wall  conditions  were  used.  The  length  of 
the  flat  plate  and  the  freestream  values  of  Mach  number, 
temperature,  Reynold's  number  and  the  turbulence  intensity  were 
(Fig.  43) ; 

L-  0.9144m,  M  -2.0,  T,  -310.93°i<-,  i^e^.  -  1 .  OK+06,  -  0.005  . 

It  was  felt  that  better  solution  comparisons  could  be  made  if 
the  same  grid  were  used  for  all  CFD  solutions  in  this  validation 
case,  so  the  rectangular  mesh  used  was  102x82  (i,j),  and  was 

provided  the  author  by  Rizzetta  [68] .  The  grid  had  a  constant 
minimum  wall-spacing  of  2.286E-05  m.  The  spacing  in  the  x- 
direction  was  a  constant  9.144E-03  m.  A  representative  grid  is 
shown  in  Fig.  44. 

4. 2. 1.1  CFD  Mach  2.0  Flat  Plate  Results.  Initial 

comparisons  of  — ,  Jf.,  and  — ,  plotted  against  the  distance  in  the 

P  u.  P. 

normal  direction,  are  shown  in  Figs.  45,  46  and  47.  The  data  for 
these  plots  were  taken  at  the  trailing  edge  of  the  flat  plate. 
Though  Fig.  45  shows  a  significant  difference  in  the  ratio  of  local 
eddy  viscosity  to  molecular  viscosity.  Figs.  46  and  47  show  that 


91 


this  difference  does  not  equate  to  as  large  a  difference  in  the 
outcome  of  the  normalized  u  and  p. 

Concern  regarding  such  a  difference  in  any  of  the  parameters 
at  the  trailing  edge  led  to  a  check  of  Cf,  the  coefficient  of 
friction.  The  plot  of  Cf  against  the  distance  along  the  flat  plate 
(Fig.  48)  showed  a  large  difference  in  the  numerical  transition 
point  (as  opposed  to  the  physical  one)  for  the  different  numerical 
models,  which  could  indeed  lead  to  the  differences  seen  above  in 
the  flow  parameters  at  the  trailing  edge. 

A  comparison  of  Fig.  49  with  Fig.  50  [47],  both  plots  of  Cf  vs 
Re^,  shows  that  the  transition  for  the  Y-S  model  occurs  within  the 
specified  transition  range  of  Re^  of  2x10^  and  6x10^  of  the  theories 
of  Prandtl  and  Blasius.  Note  that  Fig.  50  is  a  plot  for  the 
incompressible  flow  about  a  flat  plate  at  zero  degrees  incidence; 
however,  since  the  Mach  number  for  this  case  is  2.0,  there  should 
be  no  significant  change  of  Cf  due  to  compressibility,  nor  a 
significant  change  in  the  curves  in  Fig.  50  [47]  .  Once  again 
recall  that  this  is  merely  a  numerical  transition  and  not  the 
actual  physical  transition  normally  associated  with  turbulence. 

Finally,  a  plot  of  u^  vs  at  the  trailing  edge  of  the  plate 
(Fig.  51)  shows  extremely  good  comparison  between  the  two  k-e 
models  up  to  approximately  a  of  10.  From  y^  of  10  to  the  edge 
of  the  boundary  layer,  a  noticeable  difference  exists  between  the 
two  k-G  models. 
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Although  there  were  some  good  comparisons  between  the 
different  turbulence  models  and  all  of  the  models  had  been 
validated,  a  question  still  remained  as  to  whether  the  L-S  or  the 
Y-S  model  would  more  accurately  predict  the  outer  layer  of  a 
turbulent  boundary  layer.  The  logical  way  to  determine  the  k-G 
model  to  use  in  the  solver  was  to  test  them  both  on  an  experimental 
flat  plate  and  compare  the  results.  This  need  defined  the 
subsequent  step  of  the  validation. 

4.2.2  Validation  of  Y-S  k-e  Model  Against  Experimental  Data. 

Once  the  results  of  the  CFD  Mach  2.0  flat  plate  in  Section  4. 2. 1.1 
were  obtained,  there  remained  a  question  with  regard  to  the  Y-S 
model  as  to  whether  or  not  its  solution  provided  a  better 
representation  of  the  turbulent  boundary  layer  than  the  L-S  model. 
It  was  thus  decided  to  compare  the  solution  obtained  with  the  Y-S 
model  against  a  solution  obtained  with  the  L-S  model,  both  using 
the  Roe-scheme,  and  compare  both  against  well-validated  and  well- 
documented  experimental  results.  Again,  the  code  in  which  the 
turbulence  models  were  incorporated  was  perfect  gas,  full  NS,  using 
a  2-D  coordinate  system. 

The  experiment  chosen  for  validation  of  the  Y-S  k-G  model  was 
that  of  Kussoy  and  Horstman  due  to  the  fact  that  it  was  a  well- 
validated  experiment  [56] .  It  consisted  of  a  sharp  leading-edge 
flat  plate  at  -2  degrees  a  (Fig.  52)  in  an  approximately  Mach  8.2 
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flow  [56] .  The  length  of  the  flat  plate  and  the  freestream  values 
of  Mach  number,  temperature,  Reynold's  number  and  the  turbulence 
intensity  were: 

L-l.STOm,  M  -8.18,  r_.81.0“J!r,  -  7 . 938E+06,  -0.01 

The  wall  temperature  was  assumed  fixed  at  300  °K. 

The  mesh  used  was  200x90  (i,j)  and  the  minimum  wall-spacing 
was  5.8928E-06  m  (y"  =  O[l-10].  The  spacing  in  the  x-direction 
began  at  5.893E-03  m  and  increased  continually  to  1.420E-02  m.  A 
representative  grid  is  shown  in  Fig.  53. 

4.2.2. 1  Experimental  Mach  8.18  Flat  Plate  Results. 

The  plot  of  u'^  vs  y*  (Fig.  54)  did  not  initially  show 
particularly  good  agreement  with  the  computational  data  although  it 
did  exhibit  the  proper  trends.  The  lack  of  agreement  was  due  to 
the  value  of  u^.  This  value  was  originally  determined  from  the 
flow  parameters  of  the  computer  solution,  and  since  the  numerical 
transition  point  is  almost  assuredly  different  than  the 
experimental  one,  the  wall  stress  can  be  markedly  different  between 
a  computed  solution  and  experimental  data.  Once  the  value  of  u^ 
from  the  experiment  was  used  [56]  the  agreement  between  the 
computational  data  and  the  experimental  data  was  much  improved 
(Fig.  55)  .  It  can  be  seen  that  the  trends  associated  with  the 
experimental  data  are  captured  much  better  by  the  Y-S  model. 


94 


Therefore  the  Y-S  model  was  the  k-e  model  chosen  to  be  carried 


forward  for  use  in  further  turbulence  model  comparisons. 

4.2.3  Validation  of  Turbulence  Models  Against  WL  Mach  12 
Nozzle.  The  two  turbulence  model  subroutines  (B-L  and  Y-S)  were 
subsequently  incorporated  into  the  perfect  gas  and  NECR  codes  for 
validation  against  an  axisymmetric  nozzle  case  to  determine  the 
impact  of  high  temperature  effects  on  the  solution,  as  well  as  to 
determine  which  turbulence  model  provided  the  best  results.  The 
validation  case  was  the  WL  Mach  12  nozzle  experiment  [22],  with 
boundary  conditions  previously  presented  in  Sects.  3. 1.1-3. 1.4. 
Following  discussions  with  Korte  and  Scaggs  [53,63],  the  boundary 
layer  was  specified  to  be  fully  turbulent  at  the  inflow.  Should 
solutions  obtained  using  this  boundary  condition  indicate  that  it 
was  erroneous,  the  incorporation  of  some  transition  model  would  be 
required.  The  turbulence  intensity  (T^)  for  the  first  two  columns 
of  cells  was  0.01  for  the  Y-S  k-e  model. 

Both  turbulence  models  (B-L  and  Y-S)  were  run  in  each  solver. 
Both  solvers  were  run  in  perfect  gas  mode;  the  NECR  solver  was 
additionally  run  in  nonequilibri\Hn  mode. 

Two  separate  regions  were  used  to  calculate  the  entire  nozzle 
flowfield,  as  illustrated  in  Fig.  56.  The  reason  for  splitting  the 
nozzle  domain  into  two  regions  was  to  minimize  the  use  of  the 
prohibitively  small  time-marching  step  associated  with  the 
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relatively  small  subsonic  section.  The  first  region,  labelled  Item 
A,  had  a  mesh  81x126  (i,j)  and  consisted  of  the  subsonic  convergent 
section,  the  cylindrical  throat,  and  a  small  portion  of  the 
supersonic  divergent  section.  The  first  point  of  calculation  at 
the  inflow  was  located  at  an  x-location  of  9.228E-03  m  and  the 
last  point  at  the  outflow  was  located  at  an  x-location  of  3.988E-02 
m.  The  x-spacing  on  the  first  grid  varied,  but  the  nozzle  wall 
coordinates  for  the  entire  nozzle  are  presented  in  Appendix  A.  The 
second  region,  labelled  item  B,  began  10  columns  of  cells  before 
the  outflow  of  the  first  grid,  ended  shortly  before  the  nozzle 
exit,  and  was  computed  on  a  mesh  of  242x126.  For  the  second  grid, 
the  first  point  of  calculation  at  the  inflow  was  located  at  an  x- 
location  of  3.643E-02  m  and  the  last  point  at  the  outflow  was 
located  at  an  x-location  of  2.950  m.  The  value  for  y\  for  both 
grids  was  on  the  order  of  1,  achieving  a  value  of  2-3  in  the  throat 
region  of  the  first  grid  and  decreasing  after  that. 

In  order  to  maintain  second-order  accuracy  at  the  inflow  of 
the  second  mesh,  it  was  necessary  to  provide  two  initial  columns  of 
data  from  the  solution  of  the  first  grid.  The  data  for  these  two 
columns  of  cells  was  never  recomputed  in  the  solver  for  the  second 
mesh;  the  inflow  boundary  conditions  were  thus  fixed. 

4.2.3. 1  WL  Mach  12  Nozzle  Results.  When  the  computed 
static  pressure  along  the  wall  for  the  two  different  turbulence 
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models  was  compared  with  that  obtained  experimentally  (Fig.  57), 
the  B-L  turbulence  model  appeared  to  better  match  the  experimental 
data.  The  k-G  model,  on  the  other  hand,  displayed  increasing 
deviation  from  the  experimental  data  as  the  flow  traversed  the 
nozzle. 

Comparing  the  computed  Mach  number  and  Pitot  pressure  near  the 
nozzle  exit  with  experimental  data  for  the  perfect  gas  case  (Figs. 
58  and  59)  ,  the  B-L  model  once  again  provided  a  much  better 
comparison  than  the  k-G  model.  Though  the  B-L  model  appeared  to 
be  producing  a  slightly  thicker  boundary  layer  than  exhibited  by 
experiment,  the  k-G  model  produced  a  markedly  thinner  boundary 
layer  and  a  substantially  different  Mach  number  in  the  inviscid 
core.  The  Mach  number  obtained  using  the  B-L  model  was  very  close 
to  that  of  the  experiment. 

Upon  further  investigation  of  the  difference  in  the  thickness 
of  the  boundary  layer  for  the  two  turbulence  models,  an  unexpected 
change  in  the  thickness  trend  was  discovered:  the  k-G  model  was 
producing  a  thicker  boundary  layer  than  the  B-L  model  at  the  throat 
exit  location  (Figs.  60  and  61) .  The  k-G  model  thus  appeared  to 
be  more  diffusive  in  the  subsonic/sonic  section  of  the  nozzle  and 
less  diffusive  in  the  supersonic  portion  of  the  nozzle.  The  reason 
for  this  was  unclear  and  remains  to  be  investigated. 

Comparison  of  the  computed  exit  Mach  number  for  the  perfect 
gas  and  NECR  solvers  and  their  turbulence  models  showed  no 
significant  difference  between  the  perfect  gas  and  the  NECR  case 
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(Figs.  62  and  63) .  The  lack  of  difference  between  the  perfect  gas 
and  the  nonequilibrium  cases  was  presumed  to  be  due  to  the  fact 
that  the  temperature  and  pressure  in  the  flow  never  reach  levels  at 
which  nonequilibrium  and  chemical  reaction  effects  become  of 
primary  importance  [25]  . 

4.2.4  Interim  Conclusions  Drawn.  As  a  result  of  the 
validation  efforts  in  this  section,  a  number  of  significant 
conclusions  were  drawn.  One,  the  turbulence  models  had  been  well 
validated  against  both  computational  and  experimental  results. 
Two,  the  B-L  model  would  be  used  in  continuing  the  current 
research;  the  secondary  objective  of  determining  the  turbulence 
model  which  aided  in  more  accurately  representing  the  turbulence  in 
the  testbed  nozzle  flow  had  thus  been  accomplished.  Three,  there 
was  no  need  to  use  the  nonequilibrium  code  to  further  the  research 
on  the  optimization  scheme.  However,  in  future  research,  a 
nonequilibrium  form  of  the  code  will  be  developed  to  be  used  in 
high  temperature,  high  pressure  nozzles  as  the  perfect  gas  solver 
is  limited  in  such  hypersonic  nozzle  flows. 

There  were  three  additional  important  conclusions  that  were 
drawn  at  this  point.  First,  with  the  use  of  the  B-L  model  the  flow 
solution  would  provide  more  accurate  nozzle  data.  Second,  the 
assumption  that  the  flow  was  fully  turbulent  from  the  inflow  had 
been  corroborated  with  the  excellent  comparative  results  between 
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the  fully  turbulent  flow  generated  with  the  B-L  model  and  the 
experimental  data.  And  third,  the  use  of  the  B-L  turbulence  model 
would  ensure  a  more  efficient  solver,  thus  contributing  to  the 
primary  objective  of  this  research. 

4.3  Development/Validation  of  Perfect  Gas  PNS  Code  vith  Tvo 
Turbulence  Models  Incorporated 

As  mentioned  previously,  the  template  for  the  PNS  scheme  was 
that  proposed  by  Korte  [49]  due  to  the  fact  that  this  was  a  very 
thorough  source  of  PNS  formulation  using  the  Roe  scheme.  The 
upwind  flux  scheme  used  a  central  flux  approximation  for  the 
inviscid  fluxes  and  the  viscous  fluxes  were  centrally-differenced. 
Due  to  the  manner  in  which  the  computational  flux  Jacobian  was 
computed,  Vigneron's  coefficient  (co.  Eg.  49)  was  lagged  and  based 
on  the  last  complete  solution  of  the  previous  column  of  grid  cells. 

4.3.1  Validation  of  Perfect  Gas  PNS  Solver ^  Laminar  2-D  Case. 

A  suitable  validation  case  for  the  perfect  gas  PNS  solver  was 
needed,  and  Korte 's  supersonic  laminar  flow  over  a  flat  plate  was 
chosen  [49]  .  With  this  case  the  PNS  solver  could  be  validated 
against  a  well-published  test  case  and  its  solution. 

In  addition  to  validation,  the  test  case  was  run  to  determine 
if  the  GCL  terms  were  or  were  not  needed;  an  important  question 
remained  as  to  whether  a  finite  volume  scheme  required  the  extra 
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GCL  terms  which  Korte  had  outlined  [49] .  Note  that  the  GCL  terms 
are  only  necessary  if  the  grid  is  varying  in  the  y-direction  as 
space-marching  occurs.  Also  note  that  with  a  nozzle  the  y-spacing 
will  vary  in  the  axial  direction,  and  so  it  was  most  important  to 
determine  the  necessity  of  the  GCL  terms. 

The  conditions  for  the  flat  plate  validation  case,  including 
the  freestream  Mach  number  and  temperature,  the  wall  temperature, 
the  Reynold's  number  and  the  length  of  the  plate  (Fig.  64),  were 
the  following: 

M  ■  2.0,  T  -  T  -  222. 0“K,  Re^  -  1.65E*06,  L  •  1.0m. 

Two  grids  were  used  in  order  to  determine  the  necessity  of 
inclusion  of  the  GCL  terms.  The  first  was  a  rectangular  mesh  which 
varied  in  spacing  in  the  x-direction  but  not  in  the  y-direction; 
this  grid  was  used  because  it  precluded  the  need  for  the  GCL  terms . 
This  grid  is  labelled  Grid  A  and  shown  in  Fig.  65.  The  second  was 
a  grid  which  varied  in  spacing  in  both  the  x-  and  y-directions, 
being  clustered  closer  at  the  leading  edge  in  both  directions  and 
widening  in  both  directions  as  the  solution  marched  towards  the 
back  of  the  plate.  This  grid  was  used  to  determine  if  the  GCL 
terms  were  necessary,  with  such  varied  mesh  spacing,  to  generate  a 
flowfield  solution  which  compared  favorably  with  the  first  grid. 
This  second  grid  is  labelled  Grid  B,  as  shown  in  Fig.  66. 

The  spacing  in  the  x-direction  was  initially  held  to  a 
constant  value  of  l.OE-06  until  the  solution  reached  the  point 
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where  x  ^  0.025m,  at  which  point  the  spacing  in  the  x-direction  was 
allowed  to  begin  increasing  to  the  CFL  specified  space  marching 
step.  The  space-marching  step  was  initially  held  constant  to 
preclude  instabilities  associated  with  the  stagnation  region  of  the 
flat  plate  leading  edge.  However,  once  allowed  to  vary  the  spacing 
in  this  direction  was  not  allowed  to  grow  by  more  than  0.01%  at 
each  space-marching  step,  lest  instability  should  again  occur.  The 
spacing  in  the  y-direction  [49]  was  determined  from 
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The  value  of  3  was  held  constant  at  1.001  for  both  grids  in  order 
to  capture  the  boundary  layer  and  to  allow  for  a  better  comparison 
at  the  flat  plate  trailing  edge.  For  the  fixed  grid,  the  flat 
plate  was  located  at  y=0  and  the  freestream  far  edge  was  located 
at  y  =  0.8  m.  For  the  varying  grid,  the  flat  plate  was  located  at 
y  =  0  and  the  freestream  far  edge  was  located  at  y  =  0.6  m  at  the 
leading  edge  and  y  =  0.8  m  at  the  trailing  edge.  The  flow 
parameters  were  taken  from  the  x  =  0.93m  station. 
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4.3. 1.1  CFD  Laminar  Mach  2.0  Flat  Plate  Results.  Solutions 


for  Grid  A  and  Grid  B  with  and  without  GCL  terms  were  obtained. 
Comparing  the  respective  plots  of  —  vs  y  against  the  results  of 
Korte  for  the  case  excluding  GCL  terms  showed  poor  agreement 
between  the  flow  solution  obtained  using  Grid  B  and  the  other 
solutions,  particularly  in  the  peak  temperature  region  (Fig.  67) . 
With  the  GCL  terms  included,  the  plots  of  —  vs  y  compared  very 
well  for  both  grids  (Fig.  68)  .  Plots  of  —  vs  y  also  compared 
very  well,  again  with  GCL  terms  included  (Fig.  69) .  Thus  it  became 
apparent  that  the  GCL  terms  were  necessary  for  this  finite  volume 
formulation  of  PNS.  It  was  also  evident  that  the  PNS  formulation 
had  been  validated  for  the  laminar  case. 

4.3.2  Validation  of  Perfect  Gas  PNS  Solver,  Turbulent  2-D 
Case.  Next,  the  perfect  gas  PNS  solver  had  to  be  validated  using 
each  of  the  two  turbulence  models.  For  this  validation  the  Mach 
2.0  flat  plate  case  used  in  Section  4.2.1  for  the  NS  turbulence 
model  validations  was  chosen,  since  validated  solutions  and  their 
data  were  readily  available  and  required  no  additional 
computations.  The  grid  was  varied  just  as  it  was  in  the  second 
mesh  of  the  laminar  flat  plate  PNS  validation  described  in  Section 
4.3.1. 

Recall  that  the  test  case  was  a  flat  plate  at  zero  degrees  a. 
[68] .  The  length  of  the  flat  plate  and  the  freestream  values  of 
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Mach  number,  temperature,  Reynold's  number  and  the  turbulence 
intensity  were  (Fig.  43) : 

L  -  0.9144m,  M  -2.0,  r_  -  310 . -  1 . 0fJ*06,  T.  -  0.005  . 

Adiabatic  wall  conditions  were  used. 

The  validation  case  consisted  of  running  the  PNS  code  with  the 
individual  turbulence  models  (B-L  and  Y-S)  activated,  followed  by 
comparing  resultant  solutions  against  previously  obtained  NS 
results  for  the  same  turbulence  model.  The  results  would  be  used 
to  verify  that  the  PNS  code  was  computing  turbulent  flow  correctly. 

4. 3. 2.1  CFD  Turbulent  Mach  2.0  Flat  Plate  Results.  The 

first  case  run  was  that  of  the  B-L  model.  A  plot  of  normalized  u 
velocity  vs  normal  distance  from  the  plate  and  a  plot  of  Mach 
number  vs  normal  distance  from  the  plate  indicate  very  similar, 
well-matching  solutions  (Figs.  70  and  71).  It  must  be  remembered 
that  the  NS  solution  is  obtained  on  a  rather  coarse  mesh  in 
comparison  with  the  PNS  solution  in  the  streamwise  direction  (Fig. 
27),  and  that  the  fluid  center  of  a  grid  cell  is  not  necessarily 
the  geometric  center  of  the  grid  cell.  Nonetheless,  because  the 
geometric  center  is  most  easily  calculated  and  the  other  is  not, 
the  geometric  center  was  the  approximated  grid  point  associated 
with  the  flow  parameters.  This  approximation  likely  leads  to  the 
slight  disparity  between  the  two  solutions. 


103 


The  second  case  run  was  that  of  the  Y-S  model.  Again,  a  plot 
of  normalized  u  vs  normal  distance  and  a  plot  of  Mach  number  vs 
normal  distance  indicates  very  similar,  well-matching  solutions 
(Figs.  72  and  73).  The  k-e  model  appeared  to  match  even  better 
with  the  NS  solution  than  did  the  B-L  model.  This  is  possibly  due 
to  the  fact  that  the  turbulence  is  emulated  in  the  k-e  model  using 
conserved  variables,  whereas  this  is  not  the  case  with  the  B-L 
model.  Nonetheless,  the  B-L  model  still  provided  a  fairly  good 
match  between  the  two  solvers  and  it  remained  the  turbulence  model 
of  choice  in  furtherance  of  the  research. 

4.3.3  Perfect  Gas  PNS  Code  Validation  Against  Mach  12  Nozzle. 

Validation  of  the  perfect  gas  PNS  solver  on  the  WL  Mach  12  nozzle 
was  subsequently  performed.  Having  previously  validated  the 
perfect  gas  NS  solver  on  the  Mach  12  nozzle  using  each  of  the 
turbulence  models  as  described  in  Sect.  4.2.3,  the  results  of  this 
earlier  validation  would  provide  an  excellent  comparison  for  the 
perfect  gas  PNS  solver  on  the  same  case.  However,  a  couple  of 
minor  complications  presented  themselves  during  the  course  of  this 
work. 

The  inflow  to  the  PNS  solution  was  two  columns  of  cells  from 
the  first  grid  of  the  perfect  gas  NS  solver,  just  as  was  provided 
to  the  second  grid  of  the  perfect  gas  NS  solver  (Fig.  56)  . 
However,  too  large  a  space-marching  step  was  required  from  the 
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cell-center  of  the  NS  solution  to  the  cell-center  of  the  first  PNS 


column  of  grid  cells  (Fig.  74);  the  large  step  size  resulted  in 
instability,  which  caused  the  solution  to  blow  up.  This  was 
remedied  by  extrapolation  of  the  conserved  variables  using  a  finite 
difference  approximation.  The  extrapolation  was  made  using 
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The  space-marching  technique  employed  for  integration  of  the 
PNS  equations  was  initially  prohibitively  slow.  In  order  to 
accurately  capture  the  boundary  layer  in  the  throat  region  of  the 
nozzle,  it  was  necessary,  as  previously  mentioned,  to  use  a  grid 
stretching  parameter  Eq.  185)  of  1.001  in  developing  the  grid 
for  the  NS  solver  (Fig.  56)  .  And  in  order  for  the  NS  and  PNS 
solutions  to  match  up  well  at  the  inflow  to  the  latter  (Fig.  74), 
it  was  necessary  for  the  PNS  solver  to  match  the  3^  used  for  the 
grid  in  the  NS  solver.  This  value  for  was  too  small  a  value  for 
most  of  the  latter  portion  of  the  nozzle,  since  the  boundary  layer 
grows  quite  large  in  that  region;  the  restriction  on  maintained 
the  first  grid  cell  point  a  very  small  distance  off  of  the  wall 
(Fig.  75),  thus  restricting  the  space-marching  step  size  (Eq.  184). 

The  grid  adaptation  subroutine  previously  mentioned  in  Sect. 
3.3  was  utilized  to  allow  the  grid  to  vary  with  the  thickness  of 
the  boundary  layer.  This  subroutine  allows  the  user  to  specify  the 
number  of  points  in  the  laminar  boundary  layer  or  the  value  for 
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a  turbulent  boundary  layer,  while  not  allowing  the  change  from  one 
space-marching  step  to  the  next  to  be  too  large.  Utilizing  this 
subroutine  and  specifying  to  have  a  value  of  5,  the  PNS  solver 
was  able  to  compute  almost  3.0  m  of  nozzle  in  approximately  30 
minutes  of  Cray  cpu  time.  This  was  vastly  preferred  over  the  hours 
of  Cray  cpu  time  that  would  have  been  necessary  without  this 
subroutine.  In  addition,  the  value  for  3^  was  bounded.  It  was 
not  allowed  to  become  smaller  than  the  initial  value  of  1.001, 
since  the  most  restrictive  region  was  the  throat  and  the  region 
upstream  [53].  Nor  was  it  allowed  to  get  any  larger  than  1.1 
[57,71]  . 

A  specific  approach  was  necessary  to  compute  the  column  of 
grid  cells  for  the  next  space-marching  step.  Due  to  the  fact  that 
the  space-marching  step.  Ax  (Eq.  184),  and  the  wall  distance,  yv,(x) 
(Fig.  75),  were  interrelated,  coupled  with  the  fact  that 
depended  on  y„(x)  (Eq.  185),  Ax  and  Ay^^^  could  not  be  computed 
simultaneously.  The  Ax  was  computed  first;  the  Ax  was  thus  lagged 
with  the  previous  value  of  Ay„aii.  Next,  the  3^  was  calculated  as 
previously  mentioned  in  Section  3.3.  And  finally,  the  y(j)  values 
for  the  next  step  were  computed,  using  Eq.  185. 

At  this  stage,  the  perfect  gas  PNS  solver  was  ready  to  be 
validated  against  the  NS  solver,  using  the  Mach  12  nozzle  as  the 
testbed.  Concurrent  with  the  validation,  a  determination  was  to  be 
made  of  how  restrictive  to  be  with  regard  to  3  in  order  to  achieve 
a  comparable  solution  to  that  obtained  using  the  NS  solver.  Recall 
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that  the  NS  solver  used  a  mesh  with  a  B  value  of  1.001  and  the  y\ 
value  was  on  the  order  of  1  (Fig.  56)  .  Two  cases  were  run,  one 
using  a  value  of  5  and  the  other  using  a  value  of  10.  Only 
the  B-L  model  was  used  in  these  calculations. 

4. 3. 3.1  WL  Mach  12  Nozzle  Results.  Plots  of  Mach  number 
vs  radial  distance  and  Pitot  pressure  vs  radial  distance  near  the 
nozzle  exit  indicate  no  appreciable  difference  between  the  various 
solutions  (Figs.  76  and  77)  .  Though  the  value  of  plays  a  large 
part  in  the  calculated  value  of  the  shear  stress  at  the  wall  [57], 
the  plots  show  that  the  difference  in  the  output  of  the  nozzle  is 
negligible  when  is  varied  from  1,  as  in  the  PG  NS  code,  to  5 

or  10,  as  in  the  PG  PNS  code. 

4.3.4  Interim  Conclusions  Dravn.  Two  very  important 
conclusions  were  drawn  at  this  juncture  with  regard  to  the  perfect 
gas  PNS  axisymmetric  nozzle  code.  First,  based  on  the  data 
presented,  it  was  concluded  that  the  y^^^^  value  had  little  effect 
on  the  nozzle  exit  flow  as  long  as  this  value  was  between  1  and  10; 
thus,  based  on  the  recommendations  of  Korte  [9],  it  was  decided 
that  the  value  to  be  used  in  the  nozzle  optimization  scheme 

would  be  7.  Second,  the  PNS  solver  had  been  well-validated,  and 
was  ready  to  be  incorporated  into  an  efficient  optimization  scheme. 
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4.4  Objectives  Satisfied  in  Chapter  4 

Only  the  one  secondary  objective  of  the  current  research  was 
satisfied  in  this  chapter.  Recall  that  this  was  the  objective 
associated  with  choosing  a  turbulence  model  to  further  the  current 
research.  Fortunately,  the  less  computationally-intensive  B-L 
turbulence  model  proved  to  provide  the  more  accurate  testbed  nozzle 
flow  when  compared  with  experimental  data.  The  B-L  model  was  thus 
chosen  to  be  the  turbulence  model  used  in  the  remainder  of  the 
current  research.  The  other  two  objectives  remained  to  be 
satisfied. 
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5.  Development  and  Validation  of  Optimization  Scheme 

and  Code 

This  chapter  presents  the  development  and  validation  of  the 
optimization  code  in  order  to  meet  the  primary  objective  of  this 
research,  as  well  as  the  other  secondary  objective. 

In  the  development  of  the  optimization  scheme,  much  effort  was 
expended  in  order  to  achieve  a  robust  computer  code.  Methods  which 
did  not  work  well  or  required  too  much  user  interaction  were  passed 
over  in  order  to  achieve  this  goal. 

The  performance  of  a  number  of  steps  was  required  prior  to 
developing  an  automated  optimization  code.  First,  the  new  nozzle 
geometry  and  design  parameters  needed  to  be  selected.  Second,  a 
manual  optimization  search  was  required,  to  determine  if  an  optimum 
solution  to  the  hypersonic  nozzle  design  problem  could  be 
efficiently  found.  The  rationale  behind  this  was  to  minimize  the 
use  of  computer  time  in  unproductive  automatic  searches.  In  this 
manual  search,  the  response  functions  were  to  be  determined,  the 
response  surface  was  to  be  mapped,  and  the  optimization  method  was 
to  be  selected.  Following  successful  completion  of  the  manual 
optimization  search,  the  optimization  scheme  was  automated, 
validated,  and  tested  for  robustness  against  the  Mach  12  nozzle. 
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5.1  Design  Parameters  and  Nev  Nozzle  Wall  Geometry 

It  was  now  necessary  to  develop  the  design  parameters  which 
would  be  used  to  obtain  a  new  nozzle  wall  geometry.  As  previously 
mentioned,  the  two  variable  design  parameters  chosen  to  be  used  in 
the  optimization  scheme  were  the  nozzle  attachment  angle,  Sattach^ 
and  the  nozzle  exit  angle,  (Fig.  13)  .  These  two  design 
parameters  had  been  used  successfully  by  Doty  in  designing  2-D  high 
Mach  propulsive  nozzles  with  MOC/BL  techniques,  and  it  was  felt 
that  the  two  parameters  would  thus  be  sufficient  to  design  an 
axisymmetric  high  Mach  wind  tunnel  nozzle  [58] . 

A  number  of  nozzle  parameters  were  specified  and  fixed  in 
order  to  minimize  the  number  of  design  parameters  used  in  the 
optimization  procedure  (Fig.  13) .  The  x-  and  y-coordinate  values 
for  the  old  (WL  Mach  12  nozzle)  expansion  corner,  x^.^  and  as 
well  as  the  original  expansion  angle,  were  0.03279m,  0.00518m, 
and  8.842°,  respectively.  The  x-  and  y-coordinates  of  the  nozzle 
exit,  Xgxit  and  were  fixed  at  constant  values  of  3.048m  and 
0.25365m  since  these  were  the  coordinates  of  the  WL  Mach  12 
nozzle.  The  radius  of  the  attachment  corner,  was  fixed  as 
well  at  a  constant  value  of  0.00254  meters,  since  this  value 
appeared  to  well  represent  the  existing  nozzle  radius  at  the 
expansion  corner,  yexit* 

Once  ©attach  and  ©g^it  were  prescribed,  the  nozzle  wall  geometry 
could  be  obtained  using  mathematical  equations.  First,  the  x-  and 
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y-coordinates  of  the  reference  point,  and  (Fig.  13)  were 
determined  using  geometry  to  obtain 

,  (202) 

and 

Fa  =  Foe'  (203) 


where 


y  -  y  *x  .  (1-cos  (0.) )  , 

nc  oc  arc  '  ’ 

1>2  =  f -O’:- 
*1 - 2 - 


(204) 

(205) 


(206) 


(207) 


and 

dx  -  tan(02)  .  (208) 

Next,  the  X-  and  y-coordinates  of  the  attachment  point,  x^  and  y^ 
(Fig.  13),  were  determined  from  [69] 

^i,  =  ^a*-^circSin(e,,,,,„)  (209) 


and 
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=  i'a*-^clrc(l-COS(©.t,ach))  * 


(210) 


For  the  nozzle  wall  (Fig.  13)  from  the  attachment  point,  b,  to 
the  nozzle  exit,  the  wall  was  generated  using  the  equation  of  a 
skewed  parabola.  The  skewed  parabola  equation  is  given  by  [69] 

[y^-xtanT]2*qx+C2y^*C3  -  0,  (211) 


where 


(212) 

(213) 

^3  ~  ■  ^1  > 

(214) 

^2  -  -  (y«xit-^e.ittan(T))2 

(215) 

-  -(yj,--^itan(T)  )S 

(216) 

and, 

B3  =  2  [y^-x^tan{T)  ]  [tan  (T ) -tan  ]  .  (217) 

This  equation,  Eq.  211,  was  solved  using  a  Newton  iteration  in 
order  to  obtain  C^,  C2,  C3,  and  t  [44]  .  Then,  the  nozzle  wall  ra- 
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dius,  y„  (Eq.  211)  could  be  easily  generated,  given  an  axial 
location  in  the  nozzle,  x  (Fig.  13) . 

Thus,  for  the  different  solvers,  the  wall  radius,  y„,  was 
calculated  using  different  equations,  depending  on  what  the  axial 
location  was.  If  the  current  x-location  was  in  the  throat  between 
points  o  and  a,  then  the  wall  height  was  constant  at  0.0051816m. 
If  the  x-location  were  in  the  circular  throat  region,  between 
points  a  and  b,  then  the  wall  height  was  found  using  the  geometry 
of  a  circle  and  determined  by  [69] 

^  (218) 

where  0  is  determined  from  6  =  sin'^  [  (x-x  ) /r  ,  ].  Otherwise,  if  the 
x-location  were  in  the  skewed  parabola  section  of  the  nozzle 
contour  beyond  point  b,  then  Eq.  211  above  would  be  used  with  the 
appropriate  coefficients  and  skew  angle,  i,  to  generate  the  wall 
height . 

5.2  Manual  Optimization  Search 

The  first  scheme  to  be  tried  manually  provided  only  one  NS 
inflow  to  all  of  the  PNS  solutions.  The  point  in  the  nozzle  where 
the  PNS  inflow  was  initially  provided  was  at  approximately  2/3  the 
length  of  the  throat  (Fig.  78) .  The  flow  here  was  supersonic  in 
the  inviscid  core  and  it  was  felt  that  this  inflow  would  provide 
excellent  results  since  this  was  well  upstream  of  the  nozzle  wall 
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that  was  to  be  changed  (Fig.  79)  .  However,  an  oversight  occurred 
by  overlooking  the  minimum  Mach  number,  (Eq.  50)  .  Recalling 
that  the  value  for  this  parameter  is  about  1.24  for  the  stated 
value  of  the  safety  factor,  a,  the  values  for  the  Mach  number  in 
the  inviscid  core  for  this  inflow  were  approximately  1.1.  Upon 
examination  of  the  output  data  it  was  ascertained  that  Vigneron's 
coefficient,  co  (Eq.  49) ,  was  being  activated  in  the  inviscid  core 
at  the  inflow  causing  non-conservation  of  streamwise  pressure  (Fig. 
80)  . 

It  was  thus  determined  that  the  PNS  inflow  must  come  from  the 
supersonic  divergent  section  of  the  nozzle;  the  challenge 
associated  with  this  determination  was  that  the  wall  would  change 
upstream  of  the  PNS  inflow  necessitating  a  NS  solution  for  each 
change  in  the  design  parameters  (Fig.  81)  .  Since  a  NS  solution  is 
obtained  using  a  significantly  larger  amount  of  time  than  a 
comparable  PNS  solution,  potentially  many  more  hours  of  cpu  time 
might  be  necessary. 

In  order  to  minimize  the  time  associated  with  generating  a  new 
NS  solution  for  each  wall  change,  a  smaller  NS  grid,  grid  C  (Fig. 
82),  and  associated  solution  were  extracted  from  the  converged 
solution  to  the  first  NS  grid  above.  Item  A  (Fig.  56) .  In  addition 
to  being  a  smaller  grid  thus  requiring  fewer  computations  per 
iteration,  the  flow  throughout  Grid  C  was  overall  supersonic, 
allowing  for  a  larger  time  step  per  iteration  than  would  have  been 
possible  had  the  subsonic/sonic  portion  of  the  nozzle  been 
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included.  Grid  C  was  started  at  31  points  before  the  outflow  of 
the  grid  in  Item  A.  This  placed  the  first  column  of  cells  for  this 
subgrid  at  1.998E-4  m,  approximately  halfway  back  in  the  throat 
(Fig.  81).  The  outflow  column  and  all  the  other  columns  of  cells 
in  between  remained  the  same  as  in  Item  A. 

A  new  grid  which  accounted  for  the  redesigned  nozzle  contour 
due  to  a  change  in  design  parameters  was  subsequently  obtained  in 
order  to  generate  a  new  NS  solution;  from  this  new  NS  solution 
would  come  the  inflow  for  the  PNS  solver  for  each  nozzle  contour. 
Grid  C  was  used  to  generate  Grid  D  (Fig.  83),  which  had  the  same 
grid  dimensions  as  Grid  C,  31x126  (i,j). 

A  converged  NS  solution  of  Grid  D  was  then  computed.  The  NS 
solution  to  the  first  two  columns  of  cells  in  Grid  D  were  not 
recomputed;  the  solution  for  these  was  obtained  from  the  converged 
solution  to  the  grid  in  Item  A  and  was  assiamed  fixed,  presumably  a 
good  approximation  given  that  the  cells  were  far  upstream  of  the 
corner  that  would  change  with  the  design  parameters  (Fig.  83)  . 

With  Grid  D  being  used  in  the  perfect  gas  NS  solver  and  the 
converged  solution  from  Item  A  as  the  startup  solution  to  Grid  D, 
the  perfect  gas  NS  solver  did  not  appear  to  require  any  more  than 
12000  iterations  for  a  converged  solution  with  a  new  nozzle  wall 
contour  (Fig.  84)  .  Convergence  in  this  case  meant  achieving  a  PNS 
inflow  which  varied  a  negligible  amount  with  tens  of  thousands  more 
iterations  and  a  residual  which  had  decreased  2-3  orders  of  magni- 
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tude  and  was  merely  wavering  about  some  mean  value .  Each  case  was 
run  to  18000  iterations  to  ensure  convergence. 

At  this  point,  with  all  of  the  previously  performed  flow 
solver  validations,  one  could  reliably  assume  that  a  sufficiently 
accurate  NS  solution  could  be  obtained  and  provide  the  appropriate 
inflow  to  the  PNS  solver,  which  in  turn  would  produce  a 
sufficiently  accurate  nozzle  flow  solution  as  well. 

5.2.1  Determination  of  Response  Function.  The  next  challenge 
was  to  determine  the  response  function  to  be  used  in  the 
optimization.  The  first  such  function  considered  was  a  least- 
squares  residual,  specifically  the  error  in  the  Mach  number  in  the 
inviscid  high  Mach  core  at  the  nozzle  exit  (Fig.  85)  [57] .  This 

was  the  nozzle  exit  error,  defined  by 

Er  = 

luMx.i  )".  (219) 

If  necessary,  the  nozzle  centerline  error  (Fig.  86) ,  defined  by 

(220) 

would  be  added  to  Er  to  produce 

Ert-Er*Erl,  (221) 

the  total  error  function  and  another  potential  response  function. 
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In  performing  the  siammations  in  Eqs.  219  and  220,  the 
contributions  from  certain  portions  of  the  flow  might  be  preferred 
over  other  portions.  For  example,  in  the  nozzle  exit  error.  Eg. 
219,  the  errors  from  the  inviscid  core  should  contribute  heavily  to 
the  summation,  while  the  errors  from  the  boundary  layer  should  be 
deemphasized  due  to  the  likelihood  of  this  portion  of  the  flow 
dominating  the  summation.  Similarly,  in  the  nozzle  centerline 
error.  Eg.  220,  the  errors  from  the  uniform  core  near  the  nozzle 
exit  should  contribute  heavily  to  the  summation,  while  the  errors 
upstream  of  the  uniform  core  should  be  deemphasized  for  the  same 
reason  as  the  boundary  layer  in  Eg.  219. 

In  order  to  thus  accentuate  the  contributions  of  a  particular 
portion  of  the  flow  towards  a  given  response  function,  a  weighting 
factor  was  required.  Accordingly,  the  computation  of  Er  and  Erl 
would  be  modified  by  some  normalized  weighting  factor,  producing 

Er  =  Eg”  ^  "I  (222) 


and 


Erl 


\  'Umax 
^1=1 


(223) 


to  diminish  the  error  contributions  of  the  boundary  layer  to  the 
nozzle  exit  error.  Eg.  219,  and  of  the  upstream  portions  of  the 
nozzle  to  the  nozzle  centerline  error.  Eg.  220.  The  equation  for 
the  total  error  function,  Ert  (Eg.  221),  would  be  unchanged,  but 
would  now  use  the  results  from  Eqs.  222  and  223. 
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To  ensure  that  the  code  require  a  minimum  of  user-interaction 
yet  emphasize  the  preferred  portions  of  the  flow  required  the 
assumption  that  the  designer  had  no  a  priori  knowledge  of  how  large 
the  inviscid  high  Mach  core  might  be  nor  of  how  long  the  uniform 
region  might  be  in  the  latter  portion  of  the  nozzle  (Figs.  87  and 
88)  .  Hence,  the  nozzle  exit  error  summation  (Eq.  222)  was 
performed  over  the  entire  nozzle  exit.  As  will  be  discussed  next, 
the  nozzle  centerline  error  (Eq.  223)  was  performed  only  over  the 
last  half  of  the  nozzle. 

Since  it  was  felt  that  the  uniform  core  (Fig.  88)  would  be 
predominant  only  in  the  last  half  of  any  well-designed  nozzle,  the 
nozzle  centerline  summation  (Eq.  223)  was  not  performed  over  the 
first  half  of  the  nozzle.  In  addition,  with  varying  nozzle 
geometry  the  number  of  points  which  would  contribute  to  Erl  would 
vary  due  to  the  nimiber  of  space-marching  steps,  potentially  biasing 
Erl.  Specifically,  the  size  of  and  hence  the  number  of  space¬ 
marching  steps  used  in  the  PNS  solver  depends  upon  the  spacing  at 
the  wall,  as  shown  in  Sect.  3.3,  Eq.  184.  With  a  smaller 
attachment  angle,  ©attach/  the  number  of  space-marching  steps 
increases  due  to  a  smaller  Ay^a^  (Fig.  75)  .  Conversely,  with  a 
larger  ©attach  the  niomber  of  space-marching  steps  decreases  due  to  a 
larger  Ay^aii- 

Hence  it  was  decided  that  though  the  niimber  of  space¬ 
marching  steps  would  still  depend  on  the  nozzle  geometry,  the 
number  of  steps  contributing  to  the  nozzle  centerline  error 
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summation  (Eq.  223)  would  be  fixed;  this  would  remove  the 
uncertainty  associated  with  potentially  more  points  contributing  to 
the  error  summation  and  the  potential  of  biasing  the  result.  The 
last  half  of  the  nozzle  is  divided  into  50  equal  regions  (Fig.  89)  . 
Each  time  one  of  those  regions  is  entered  with  a  space-marching 
step,  the  error  associated  with  the  current  point  on  the  centerline 
contributes  to  the  nozzle  centerline  error.  When  the  end  of  the 
nozzle  is  reached,  that  point  contributes  to  Eq.  223  as  well. 
Accordingly,  iimax  (Eq.  223)  has  a  constant  value  of  51,  which  is 
specified  in  the  code. 

Another  issue  which  arose  dealt  with  the  aforementioned  weight 
factors.  The  criteria  to  be  met  by  a  weighting  factor  were  as 
follows:  the  errors  associated  with  the  boundary  layer  at  the 
nozzle  exit  (Eq.  222)  and  the  non-uniform  portion  of  the  centerline 
(Eq.  223)  be  diminished;  the  weighting  factors'  sum  over  the 
domain  of  interest  be  equal  to  1;  the  weighting  factors  always  be 
positive;  and  the  weighting  factors  maintain  the  nozzle 
optimization  as  a  search  for  a  minimiam. 

It  was  decided  that  the  nozzle  exit  and  nozzle  centerline  x- 
coordinate  Mach  profiles  most  resembled  a  step-function  with  a 
finite  rise  time,  which  is  represented  by 

f{fc)  =  f{=0)  d-e'”)  ,  (224) 

where  t  is  a  constant  (Fig.  90)  [70].  Thus,  when  t  is  near  0,  f(t) 
is  near  zero.  And  when  t  goes  to  infinity,  f(t)  approaches  f(“). 
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An  analogous  case  exists  with 


Q  -  6Xp  (~l^ifl|ax,  j'-^deaignl  ^ 

used  in  Er,  Eq.  222,  and 


(225) 


(226) 


used  in  Erl,  Eq.  223.  In  this  case,  when  the  power  of  the 
exponential  is  near  zero,  that  is  the  Mach  number  is  near  M^gsign, 
the  exponential  term  is  close  to  1;  this  would  be  the  case  in  the 
inviscid  high  Mach  uniform  core  at  the  exit  if  the  nozzle  were 
designed  correctly  (Fig.  85) .  When  the  power  of  the  exponential  is 
not  near  zero,  that  is  the  Mach  number  is  not  very  close  to  M^esign^ 
the  exponential  term  is  close  to  zero;  this  would  be  the  case  in 
the  boundary  layer  at  the  nozzle  exit  or  the  non-uniform  portion 
along  the  nozzle  centerline  (Fig.  86)  .  Hence,  the  portions  of  the 
flow  which  are  most  important  are  emphasized  and  the  portions  which 
are  not  important  are  deemphasized.  Note  that  all  of  the 
aforementioned  criteria  for  weighting  factors  were  met  with  Eqs. 
225  and  226. 

Once  again,  these  weighting  factors  preclude  the  need  to 
specify  some  portion  of  the  solution  to  sample,  which  would 
definitely  change  with  different  nozzles.  This  contributes  to  the 
minimum  interaction  nature  of  the  computer  code. 
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The  nozzle  exit  error,  Eq.  219,  was  the  sole  response  function 
initially  used,  but  merely  having  uniform  flow  at  the  outflow  of 
the  nozzle  is  no  guarantee  that  there  is  indeed  a  uniform  inviscid 
core  with  a  minimum  of  disturbances  (Fig.  91) .  Thus  it  was  decided 
that  the  total  response  function,  Ert  (Eq.  221),  would  provide  the 
response  values,  z^  (Eq.  114),  for  the  optimization  scheme  since  it 
sought  to  minimize  both  of  the  error  functions  simultaneously. 
Additionally,  Eqs.  222  and  223  would  provide  the  input  to  Ert,  thus 
providing  a  more  user-friendly  optimization  procedure. 

5.2.2  Mapping  the  Response  Surface.  In  order  to  provide  some 
information  on  the  response  surface's  appearance,  the  design 
parameters  (©attach/  ^exit)  were  varied  systematically.  The  mapping 
of  the  response  functions  was  generated  by  specifying  5  fixed 
values  of  ©attach/  arid  each  of  these  values,  6  fixed  values  of 
©exit  (Fig.  92)  .  The  values  for  the  design  parameters  used  in  this 
study  are  presented  in  Table  1.  Each  point  in  Table  1  was  input  to 
the  NS/PNS  solver  to  obtain  the  associated  error  functions;  these 
were  then  used  to  generate  maps  of  the  individual  response  surfaces 
generated  by  using  Er  (Eq.  222),  Erl  (Eq.  223),  and  Ert  (Eq.  221). 
The  conditions  for  this  study  were  previously  presented  in  Sects. 
3. 1.1-3. 1.4  and  Sect.  3.3.1. 
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Table  1.  Sampling  matrix  for  response  surface  mapping. 


e 


attach 


0 


exit 


8.74° 

0.00° 

0.94° 

1.89° 

O 

00 

00 

• 

CM 

3.78° 

4.58° 

10.6° 

0.00° 

0.94° 

1.89° 

2.83° 

3.78° 

4.58° 

12.77° 

0.00° 

0.94° 

1.89° 

2.83° 

3.78° 

4.58° 

14.9° 

0.00° 

0.94° 

1.89° 

2.83° 

3.78° 

4.58° 

16.8° 

0.00° 

0.94° 

1.89° 

2.83° 

3.78° 

4.58° 

Geometrical  constraints  placed  physical  limits  on  the  design 
variables,  ©attach  and  ©exit*  The  lower  limit  of  ©attach  was  determined 
by  the  constraint  of  the  cone  angle,  seen  if  a  line  were  drawn  from 
the  corner  of  the  throat  exit  to  the  wall  at  the  outflow  of  the 
nozzle  (Fig  93)  .  The  upper  limit  for  ©attach  was  determined  by  the 
constraint  of  the  angle  at  which  flow  separation  would  begin  to 
occur,  since  this  would  be  unacceptable  as  an  inflow  to  the  PNS 
solver.  The  lower  limit  of  ©^xit  was  determined  by  a  parallel  wall, 
equating  to  a  zero  degree  angle  for  the  wall  at  the  nozzle  exit 
(Fig.  94)  .  And  finally,  the  upper  limit  of  ©g^it  was  determined  by 
the  constraint  of  the  cone  angle  already  mentioned.  It  was  hoped 
that  obtaining  the  response  surface  map  would  accomplish  several 
objectives.  First,  it  would  allow  one  to  see  if  the  optimization 
scheme  was  marching  towards  a  minimum  or  not.  Second,  it  would 
provide  more  information  regarding  the  response  surface,  including 
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the  overall  appearance.  Third,  it  would  aid  in  the  validation  of 
the  optimization  procedure. 

5.2.2. 1  Results  of  Response  Surface  Mapping.  With  the 
design  parameter  coordinates  provided  in  Table  1  above,  response 
surfaces  were  generated  using  Tecplot  to  obtain  maps  for  Er,  Erl, 
and  Ert  (Figs.  95,  96  and  97).  Analysis  of  the  response  surfaces 
provided  a  number  of  insights.  First,  though  they  each  exhibited 
a  point  of  global  minimum,  the  areas  associated  with  the  global 
minima  for  Er  and  Erl  were  substantially  larger  than  that 
associated  with  Ert.  Second,  they  all  exhibited  local  minima 
different  from  the  global  minima  for  each  surface.  Thus,  the 
potential  existed  for  the  optimization  scheme  to  find  a  local 
minimum  and  not  the  global  one.  Third,  the  maps  did  not  exhibit 
simple  first-  or  second-order  behavior,  particularly  in  the  low 
^attach  region  or  the  high  ©exit  region;  indeed,  many  forms  of  behavior 
were  evident  on  the  maps,  including  ridges  and  plateaus.  It  should 
be  pointed  out  that  the  potential  exists  for  other  local  minima  not 
displayed  by  Tecplot  to  exist  due  to  the  methodology  used  within 
Tecplot  to  generate  the  map  surface. 

Two  conclusions  were  drawn  from  the  above  observations.  First, 
it  appeared  that  too  large  an  area  was  associated  with  the  global 
minima  for  Er  and  Erl .  Since  its  map  appeared  to  indicate  a  more 
precise  global  minimum  region,  this  confirmed  that  Ert  should  be 
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used  as  the  response  function.  Second,  since  the  response  surface 
map  indicated  curvature  in  the  region  of  the  local  minimum,  it  was 
decided  to  use  a  second-order  search  exclusively  in  trying  to 
locate  the  global  minimum.  Another  factor  in  this  decision  was 
that  fewer  sampling  points  were  expected  to  be  required  than  would 
be  in  the  steepest  descent  method.  It  was  expected  that  the 
second-order  search  method  would  find  the  global  minimum,  despite 
the  variety  of  surface  behavior  displayed  on  the  map. 

5.2.3  First  Scheme,  Second-Order  Response  Surface  Methodology 
(RSM)  Search.  With  the  response  functions  and  the  response  surface 
mapping  addressed  above,  the  second-order  RSM  optimization 
procedure  was  implemented.  The  starting  point  for  this  case  was 
(©attach^  ©exit)  of  (12.8°,  2.30°)  with  de^ttach  of  5°  and  dOexit  of  1.5°. 
Each  time  a  new  minimum  was  found  using  the  procedures  provided  in 
Sect.  2. 8. 2. 2,  the  new  center  was  shifted  to  that  point,  and  the 
values  of  dSattach  dSexit  were  cut  by  half  (Fig.  19)  .  The  search 
was  continued  with  successive  iterations  until  either  a  minimum  had 
been  located  or  it  was  determined  that  the  optimization  procedure 
was  not  efficiently  seeking  out  a  minimum.  The  conditions  for  this 
study  were  previously  presented  in  Sects.  3. 1.1-3. 1.4  and  Sect. 
3.3.1. 
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5.2.3. 1  First  Scheme  Results.  The  inability  of  the  second- 
order  search  method  to  locate  the  global  minimum  in  an  expedient 
manner  is  evidenced  by  the  trail  on  the  response  surface  map  (Fig. 
98)  .  After  five  iterations,  the  second-order  search  method  had 
still  not  located  the  global  minimum  and  appeared  to  be  having 
considerable  difficulty  in  doing  so.  This  confirmed  that,  with  the 
response  surface  contours,  an  exclusive  second-order  search  method 
was  not  the  best  choice  and  another  method  needed  to  be  explored. 

5.2.4  Second  Scheme,  Korte  Least  Squares  (LS)  Method.  To 

explore  potential  improvement,  the  LS  method  of  Korte  was  tried 
next,  using  the  procedures  outlined  in  Sect.  2.8.1.  The  method  of 
Korte  was  used  because  that  method  was  well  proven  and  required 
only  one  NS  solution  per  sampling  grid;  this  was  due  to  the  use  of 
extremely  small  d6’s  in  obtaining  the  sampling  grid,  as  will  be 
shown.  In  addition,  only  three  sampling  points  were  necessary  per 
iteration,  as  compared  with  at  least  nine  for  the  RSM  methods. 

To  start,  the  last  identified  minimum  of  the  RSM  second-order 
search  method  (Fig.  98)  provided  the  starting  point  (point  1)  of  a 
three  point  stencil  (Fig.  99) .  The  three  points  of  the  stencil 
always  consisted  of  the  last  minimum  found  (point  1),  and  two  other 
points  (points  2  and  3),  which  were  generated  by  alternately 
holding  one  of  the  design  parameters  of  the  starting  point 
constant,  while  the  other  was  varied  by  some  small  d9.  As 
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reconimended  in  [57],  the  original  d0  was  0.001  9.  Thus,  the  three 
points  of  the  first  stencil  consisted  of  (©attach/  ©exit)  values  of 
(14.45°,  1.800°)  for  point  1,  (14.46°,  1.800°)  for  point  2,  and 
(14.45°,  1.798°)  for  point  3.  The  value  of  d6  was  subsequently 
decreased  to  0.0001  0  and  lower  values  in  an  attempt  to  locate  a 
minimum  in  an  efficient  manner.  The  conditions  for  this  study  were 
previously  presented  in  Sects.  3. 1.1-3. 1.4  and  Sect.  3.3.1. 

The  components  of  the  Jacobian  matrix,  Eq.  99,  were  determined 
using  a  first  order  approximation,  and  the  error  functions  in  Eqs. 
222  and  223  were  used  not  as  a  sum  as  in  Eq.  221,  but  independently 
in  this  case.  Thus,  given  that  p^  is  Er,  is  Erl,  a^  is  ©attach/  and 
a2  is  ©exit/  then  the  (i,j)  component  of  the  Jacobian  matrix  would 
be  determined  by 


da.  a.  -a 

J  J«  Jt 


(227) 


where  m  is  the  index  representing  a  point  of  the  stencil  (point  2 
or  point  3,  Fig.  99) ;  m  is  2  if  j  is  1,  accounting  for  a  change  in 
©attach/  or  m  is  3  if  j  is  2,  accounting  for  a  change  in  ©exit* 

The  solution  (©attach/  ©exit)  fon  the  newly  computed  minimum  point 
could  be  obtained  once  the  Jacobian  matrix,  Eq.  99,  had  been- fully 
computed.  First,  Eq.  98  was  solved  for  using  a  simple  2x2 
matrix  inversion  procedure.  Then,  Aa^^  was  added  to  Eq.  100,  to 
obtain  the  newly  determined  minimum  point,  a^^^^ . 
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5. 2. 4.1  Second  Scheme  Results.  Though  the  Korte  LS  method 
did  find  the  global  minimum  at  one  point,  it  required  much 
interaction  with  the  code  in  order  to  maintain  d0  in  a  range  that 
would  not  cause  the  optimization  procedure  to  wander  wildly  about 
the  surface  in  search  of  a  minimum.  The  path  taken  by  the 
optimization  procedure  is  shown  on  the  response  surface  map  (Fig. 
100)  . 

Unfortunately,  local  minima  consistently  began  appearing 
between  two  very  close  sampling  points,  where  presumably  no  local 
minima  should  be  (Fig.  101) .  It  became  clear  that,  with  such  small 
dG's,  any  changes  that  occurred  in  the  response  function  were  below 
the  numerical  noise  level,  meaning  that  the  response  function 
output  was  not  to  be  trusted  to  determine  a  gradient  direction.  An 
attempt  to  raise  the  value  of  the  d0’s  led  to  erratic  search 
behavior  (Fig.  102) .  Hence,  the  attempt  to  use  the  Korte  LS  method 
was  abandoned. 

5.2.5  Third  Scheme ,  Steepest  Descent  RSM  Method.  The  next 
scheme  developed  was  the  steepest  descent  method.  The  procedures 
used  were  previously  given  in  detail  in  Sect.  2. 8. 2.1.  The 
conditions  for  this  study  were  previously  presented  in  Sects. 
3. 1.1-3. 1.4  and  Sect.  3.3.1. 

Beginning  with  the  same  center  and  d0's  used  in  the  beginning 
of  the  second-order  search  above,  the  five  point  stencil  was  used 
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to  generate  the  presumed  direction  of  steepest  descent  (Fig.  15) . 
Initial  g  values  (Eg.  121)  of  1,  ^|2,  sfZ,  2  were  tried.  The  g  values 
next  needed  to  be  decoded  into  the  (©attach/  ©exit)  coordinates  using 
Eqs.  125  and  126.  Each  g  value  was  then  run  through  the  NS/PNS 
solvers  in  order  to  determine  the  response  function  associated  with 
each.  If  the  range  of  the  initial  g  values  was  not  sufficient  to 
locate  a  response  surface  inflection  point  along  the  g  vector, 
larger  values  of  g  were  used  successively  until  an  inflection  point 
was  located  (Fig.  103)  .  For  instance,  if  the  error  functions 
associated  with  the  previously  used  g  values  indicated  no 
inflection  point,  some  amount  6g  was  added  to  the  last  used  g  value 
to  obtain  the  next  g  value  to  be  used  in  the  search.  This 
continued  until  an  inflection  point  was  bracketed  by  assigned  g 
values  (Fig.  103),  and  then  the  three  closest  g  value  error 
functions  were  used  to  obtain  Eq.  133,  and  its  (©attach/  ©exit) 
coordinates.  The  value  of  was  then  run  through  the  NS/PNS 
solvers  to  verify  that  a  new  minimum  did  exist  at  that  (©attach/  ^exit) 
coordinate.  With  the  new  minimum  found,  that  (©attach/  ©exit)  point 
became  the  center  of  the  five  point  stencil  (Fig.  11)  and  the  d©'s 
were  decreased  in  the  hopes  of  zeroing  in  on  the  minimum.  This 
procedure  was  iterated  until  a  minimum  had  been  found  which  could 
not  be  improved  upon. 
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5. 2. 5.1  Third  Scheme  Results.  Surprisingly  enough,  the 
steepest  descent  method  sought  out  the  global  minimum,  as  evidenced 
by  the  path  followed  on  the  Ert  response  surface  map  (Fig.  104) 
This  minimum  was  located  within  three  iterations,  and  two 
subsequent  iterations  failed  to  find  another  in  the  local  vicinity. 
Plots  of  the  Mach  number  vs  radial  distance  at  the  exit  and  Mach 
number  vs  axial  distance  along  the  centerline  indicate  that  the 
minimum  does  correspond  to  a  uniform  Mach  12  inviscid  core  at  the 
exit,  at  least  as  is  likely  obtainable  using  only  the  two  design 
parameters  (Figs.  105  and  106). 

5.2.6  Comparison  vith  Original  Mach  12  Nozzle  Design.  A 

comparison  between  the  solutions  for  the  new  and  old  Mach  12 
nozzles  indicates  that  the  new  nozzle  design  produces  a  more 
accurate  test  section  flow,  but  with  perhaps  more  disturbances  in 
the  flow.  Examination  of  Fig.  107  shows  the  exit  profile  is  much 
nearer  to  the  design  mach  number  for  the  new  nozzle  than  it  is  for 
the  old  nozzle  design.  Examination  of  Fig.  108  however  indicates 
more  disturbances  in  the  flow  of  the  new  nozzle  design  as  compared 
with  the  original.  This  may  indicate  that  there  are  limitations  to 
the  optimization  scheme  with  the  selected  design  parameters.  The 
optimization  scheme  does  nevertheless  optimize  the  nozzle  to  the 
best  of  its  abilities. 


129 


5.3  Automating  the  Search 

With  the  steepest  descent  method's  successful  location  of  the 
global  minimum,  the  entire  optimization  scheme  was  coded,  including 
the  NS  solver  of  the  subgrid  (Grid  D,  Fig.  83),  the  PNS  solver,  the 
steepest  descent  method,  the  second-order  search  method,  and  all  of 
the  logic  previously  mentioned.  The  conditions  for  this  study  were 
previously  presented  in  Sects.  3. 1.1-3. 1.4  and  Sect.  3.3.1.  The 
starting  point  and  the  d0's  remained  unchanged  from  the  initial 
starting  values  for  the  manual  steepest  descent  method.  However, 
the  initial  values  for  g  in  the  direction  of  steepest  descent  were  ^ ; 
the  reason  for  the  difference  in  the  initial  values  of  g,  as 
compared  with  the  manual  steepest  descent  search,  was  that  it  was 
feared  that  with  such  large  initial  values  the  automated  search 
might  immediately  violate  the  coordinate  limits  on  the  response 
surface  (Fig.  109) . 

If  a  new  minimum,  as  compared  with  the  center  point,  had  not 
been  found  among  those  three  points,  the  response  surface  was 
apparently  increasing  in  the  calculated  direction  of  steepest 
descent  (Fig.  110) .  The  optimization  scheme  then  reverted  to  the 
second-order  search,  since  5  of  the  9  sampling  point  response 
functions  had  already  been  computed  (Fig.  12) . 

If  the  response  surface  was  decreasing  in  the  direction  of 
steepest  descent,  the  optimization  scheme  then  checked  to  see  if  an 
inflection  point  had  been  bracketed  by  the  g  values  used  (Fig. 
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103)  .  If  an  inflection  point  had  been  bracketed,  the  optimizer 
then  computed  the  g  value  for  the  calculated  minimum,  (Eq. 
140)  ,  using  the  second-order  bivariate  regression  equation,  Eq. 
136,  and  computed  the  error  function  for  that  point using  the 
NS/PNS  solvers.  The  error  function  values  for  the  bracketing  g 
values  and  were  then  compared,  with  the  (©attach^  ©exit)  point 
associated  with  the  new  found  minimiim  then  becoming  the  center  for 
the  next  sampling  grid  (Fig.  15) . 

If  an  inflection  point  had  not  been  found,  the  optimizer 
continued  to  search  in  the  direction  of  steepest  descent  until  one 
had  been  located.  After  the  g  value  of  /2  had  been  passed,  the 
delta  to  the  next  g  value  from  the  current  one  was  obtained  by 
multiplying  the  previous  delta  by  1.1  (Fig.  103). 

Assuming  that  a  new  minimum  had  been  found,  the  procedure  used 
in  the  steepest  descent  method  then  decreased  the  d9's,  and  with 
the  new  minimum  as  the  center  of  the  new  five  point  stencil  began 
searching  anew  for  the  next  minimum  (Fig.  15) .  If  this  were  the 
second  iteration,  then  dGattach  would  be  decreased  to  l/5th  its 
original  value  and  dSe^it  would  be  decreased  by  1/2;  the  decrementing 
of  the  d0's  was  accomplished  in  this  manner  in  order  to  diminish 
the  height  to  width  ratio  of  the  sampling  stencil  in-  non- 
transformed  space  (Fig.  Ill) .  Additionally,  the  large  initial 
dSattach  was  due  to  the  fact  that  there  was  initially  a  much  larger 
area  to  sample  in  that  direction  (Fig.  109) .  In  subsequent 
iterations,  both  d0’s  were  decreased  to  1/2  of  their  last  values 
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when  a  new  iteration  began;  the  d6's  were  however  not  allowed  to 
decrease  to  less  than  0.1  because  it  was  felt  that  any  changes 
below  that  level  would  produce  changes  in  the  error  function  which 
would  be  attributed  to  the  numerical  noise  rather  than  to  the 
design  parameters. 

The  iterative  process  continued  until  a  minimum  had  been 
located  and  one  of  two  events  occurred.  Either  the  user  perceived 
no  discernible  difference  in  the  output  from  the  nozzle  (Mach 
number  and  Pitot  pressure  profiles)  with  the  detection  of  a  new 
minimum  or  no  other  minima  had  been  found  in  the  area  and  the  limit 
on  the  d6's  had  been  reached. 

5.4  Efficiency  of  Scheme  Methodology  Tested 

To  determine  if  a  minimiua  could  be  obtained  in  a  more 
efficient  manner,  a  viscous  Mach  12  hypersonic  nozzle  optimization 
search  was  performed  with  a  slightly  different  sampling  stencil. 
Except  as  mentioned  here,  all  of  the  nozzle  details  remain 
unchanged  from  Sect.  5.3.  The  same  initial  stencil  was  used  as 
that  of  the  search  performed  in  Sect.  5.3.  However,  the  second  and 
later  stencils  in  non- transformed  space  each  constituted  a  square, 
rather  than  a  rectangle,  with  the  newly  found  minimum  at  the  center 
(Fig.  112)  .  Thus  for  the  second  iteration,  d0 attach  was  decreased  to 
1/5  of  its  value  and  dSe^it  to  2/3  of  its  value  in  order  to  square 
up  the  sampling  stencil.  In  subsequent  iterations,  the  dO's  were 
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decreased  to  2/3  of  their  previous  values,  rather  than  1/2,  as  in 
Sect.  5.3.  With  these  exceptions,  the  search  was  performed  as  it 
was  described  in  Sect.  5.3. 

5.5  Optimization  Code  Tested  for  Robustness 

In  addition  to  verifying  that  the  optimization  code  had  been 
developed  to  a  satisfactory  level  of  efficiency,  it  was  necessary 
to  ensure  that  it  was  also  at  a  satisfactory  level  of  robustness. 
In  other  words  it  was  important  to  ensure  that  given  a  number  of 
completely  different  starting  points  and  non-transformed  stencils, 
the  same  global  minimum  could  be  obtained.  This  would  validate 
that  the  finding  of  the  global  minimum  had  not  been  case-specific. 

Two  separate  test  cases  composed  of  two  different  starting 
points  were  run.  All  of  the  conditions  for  the  following  tests  are 
the  same  as  those  described  in  Sect.  5.3. 

Note  that  the  seemingly  less  efficient  sampling  grid  reduction 
was  used  in  the  following  two  test  cases.  This  was  due  to  the  fact 
that  the  computations  for  these  cases  were  performed  simultaneously 
as  an  efficiency  check  on  sampling  grid,  and  it  was  expected  that 
the  square  non-transfomed  sampling  grid  and  lesser  decrementing  of 
the  d0's  would  be  more  efficient.  As  seen  above,  this  did  not 
appear  to  be  the  case.  Although  the  change  in  sampling  grid 
reduction  resulted  in  more  run  time,  the  objective  of  the  tests 
were  not  compromised. 
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5.5.1  First  Test.  The  first  test  objective  was  to  determine 
if  the  global  minimum  could  be  obtained  with  a  different  starting 
point.  The  search  was  begun  with  a  starting  point  of  (©attach^  ©exit) 
of  (15,  1),  which  was  towards  the  maximum  for  ©attach  arid  the  minimum 
for  ©exit  (Fig.  109).  The  d©'s  were  (2.5,  1),  and  the  same  dec¬ 
rementing  occurred  as  in  the  previous  search,  where  d©attach  '^as 
decreased  to  1/5  of  its  value  and  d©exit  to  2/3  of  its  value  for  the 
second  iteration.  For  subsequent  iterations,  the  decrements  in  the 
d©'s  were  to  2/3  of  their  previous  values. 

5.5.2  Second  Test.  The  second  test  for  the  optimizer  had  the 
objective  of  determining  if  the  global  minimum  could  be  obtained 
starting  from  a  point  which  was  close  to  a  local  minimum  but  not 
the  global  minimum.  The  starting  point  for  this  search  was  located 
at  (©attach^  ©exit)  of  (10,  3.5),  which  was  located  near  the  minimum 
for  ©attach  arid  the  maximum  for  ©exit*  seen  in  Fig.  97,  this  point 
is  closer  to  the  local  minimum  located  at  approximate  values  of 
(©attach/  ©exit)  of  (12.5,  4.0).  The  d©'s  had  values  of  (1.0,  1.0). 
The  decrements  in  the  d©'s  were  to  2/3  of  their  previous  values  for 
all  iterations. 

5. 6  Use  of  NECR  Solver  to  Justify  Use  of  Perfect  Gas  Solver 

It  was  still  necessary  to  prove  that  for  these  operating 
conditions  a  hypersonic  nozzle  of  this  type  would  not  require  the 
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use  of  a  nonequilibrium  solver  in  order  to  correctly  optimize  the 
nozzle  contour.  Thus,  once  the  newly  optimized  nozzle  shape  was 
obtained,  it  was  necessary  to  obtain  a  solution  of  the  new  nozzle 
with  the  aid  of  the  nonequilibrium  solver.  All  other  conditions 
are  as  previously  specified  with  the  exception  of  the  wall  contour 
of  the  testbed  WL  Mach  12  nozzle. 

5.  7  Off-Design  Conditions 

The  computational  capability  of  testing  off-design  conditions 
is  expected  to  be  of  interest  in  nozzle  optimization.  Therefore, 
as  a  further  test  of  the  perfect  gas  NS  code  two  off-design 
conditions  were  run.  The  first  case  was  that  of  increased  exit 
pressure.  The  exit  pressure  at  every  point  along  the  nozzle  exit 
was  set  at  twice  that  of  the  nozzle  static  pressure  at  the 
centerline  when  no  back  pressure  was  specified,  as  was  previously 
the  case  in  the  NS  solution  for  the  Mach  12  nozzle  (Sect.  3.1.4). 
This  situation  might  be  analogous  to  a  buildup  of  pressure  in  the 
vacuum  reservoir  of  a  wind  tunnel  or  to  a  slight  blockage  of  the 
flow  due  to  the  model.  The  second  case  was  that  of  decreased  exit 
pressure.  The  exit  pressure  was  set  to  half  that  of  the  nozzle  at 
the  centerline  when  no  back  pressure  was  specified.  This  situation 
might  simulate  a  startup  of  a  wind  tunnel.  Note  that  these  cases 
could  not  be  run  using  the  perfect  gas  PNS  code,  since  there  is  no 
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capability  of  specifying  a  back  pressure  using  PNS  in  the 
methodology  utilized  in  this  research. 
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6.  Results  and  Conclusions  from  Validation  of 


Optimization  Code 

This  chapter  presents  the  results  of  the  automated  validation 
cases  and  the  off-design  cases  from  Ch.  5  (Sects.  5. 3-5. 7)  and 
addresses  conclusions  drawn  and  objectives  met  or  not  met.  Recall 
that  the  primary  objective  of  this  research  was  to  develop  and 
validate  a  computer  code  which,  with  Navier-Stokes  (NS)  and 
Parabolized  Navier-Stokes  (PNS)  solvers  incorporated,  would 
accurately  and  robustly  redesign  the  testbed  nozzle  utilizing  a 
simple  optimization  scheme,  given  the  fixed  length,  throat  radius 
and  exit  radius  of  the  original  Mach  12  nozzle. 

To  this  point,  the  primary  objective  and  the  secondary 
objective  of  running  the  off-design  cases  on  the  newly  designed 
nozzle  had  yet  to  be  achieved.  However,  the  NS  and  PNS  flow  codes 
had  been  well-validated.  In  addition,  the  Baldwin-Lomax  (B-L) 
model  had  proven  to  be  the  turbulence  model  of  choice,  satisfying 
one  secondary  objective.  Also,  the  flow  codes  had  been  used  in  the 
manual  operation  of  the  optimization  scheme  and  the  global  minimum 
had  been  achieved  (Sect.  5.2)  .  Finally,  the  stage  was  set  for  the 
testing  of  the  automatic  operation  of  the  optimization  code  with 
the  flow  solvers  incorporated. 
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6.1  Results  of  Automating  the  Search 

The  optimization  scheme  once  again  found  the  global  minimum, 
and  in  this  case  successful  automation  required  four  to  five 
iterations  to  achieve  the  global  minimum.  By  four  iterations,  the 
response  function  minimum  had  attained  a  value  of  0.02986  as 
compared  with  0.1305  for  the  first  iteration.  The  fifth  iteration 
delivered  a  value  for  the  minimxim  of  0.02946.  The  search  was  then 
terminated  as  the  d0's  had  reached  their  limit  of  0.1.  A  plot  of 
the  path  taken  by  the  optimization  scheme  is  imposed  on  the 
response  surface  map  in  Fig.  113. 

In  order  to  observe  the  difference  associated  with  each  of  the 
three  response  function  minima  discussed  above,  plots  of  Mach 
number  vs  radial  distance  at  the  nozzle  exit  and  Mach  number  vs 
axial  distance  along  the  nozzle  centerline  are  shown  (Figs.  114  and 
115).  While  there  was  a  noticeable  difference  between  the  plot  for 
the  first  minimum  and  the  other  two,  there  is  no  discernible 
difference  between  the  plots  for  the  last  two  minima  found. 

Recall  that  though  the  initial  sampling  stencil  was  the  same 
as  used  in  the  manual  steepest  descent  search,  the  stencil  along 
the  path  of  steepest  descent  was  made  different  to  avoid  the 
potential  of  violating  one  of  the  design  parameter  boundaries  (Fig. 
109)  .  This  conservative  approach  appears  to  have  added  one  to  two 
iterations  in  the  search  to  achieve  the  global  minimum.  Complete 
resultant  path  data  is  presented  in  Appendix  B. 
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6.2  Results  of  Efficiency  Test 

As  suspected,  the  sampling  stencil  did  have  a  minor  effect  in 
the  efficiency  of  the  optimization  scheme.  The  new  sampling  method 
led  to  an  increase  in  the  number  of  iterations  necessary  to  attain 
the  global  minimum.  This  case  required  six  iterations  to  reach  the 
global  minimum,  as  opposed  to  the  four  to  five  iterations  seen 
earlier.  The  path  taken  by  the  optimization  scheme  is  shown  on  the 
response  surface  map  (Fig.  116) .  The  original  non-transformed 
rectangular  stencil  appears  to  provide  quicker  convergence,  as  did 
the  more  aggressive  decrementing  of  the  d6*s.  Complete  resultant 
path  data  is  presented  in  Appendix  B. 

6.3  Results  of  First  Test  for  Robustness 

The  global  minimum  was  attained  using  this  new  starting  point. 
Five  iterations  were  required  to  attain  the  global  minimum,  and  the 
search  required  the  code  to  revert  to  the  second-order  search 
twice.  The  path  taken  by  the  optimization  scheme  is  shown  plotted 
on  the  response  surface  map  (Fig.  117) .  Complete  resultant  path 
data  is  presented  in  Appendix  B. 

6.4  Results  of  Second  Test  for  Robustness 

The  second  new  starting  point  did  not  locate  the  global 
minimum.  It  did  locate  a  local  minimiam,  but  not  the  local  minimum 
towards  which  it  appeared  likely  to  be  attracted.  The  local 
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minimum  found  did  not  appear  on  the  response  surface  map, 
presumably  due  to  the  sparseness  of  the  mapping  grid  and/or  the 
methodology  used  by  the  Tecplot  software  package  to  interpolate 
between  points.  This  case  required  3  iterations  to  attain  the 
local  minimum,  reverting  to  second-order  search  twice.  The  path 
taken  by  the  optimization  scheme  is  shown  plotted  on  the  response 
surface  map  (Fig.  118) .  Complete  resultant  path  data  is  presented 
in  Appendix  B. 

6.5  Results  of  NECR  Solver  vith  Netr  Nozzle  Design 

The  new  nozzle  design  produced  no  significant  differences 
between  the  perfect  gas  and  the  NECR  NS  solvers .  Plots  taken  from 
near  the  nozzle  exit  of  Mach  number  vs.  radial  distance  and  of 
Pitot  pressure  vs.  radial  distance  are  presented  in  Figs.  119  and 
120.  A  plot  of  Mach  number  vs.  axial  distance  along  the  nozzle 
centerline  is  presented  in  Fig.  121.  Based  on  these  results,  it  is 
concluded  that  the  use  of  the  NECR  solver  in  optimizing  the  testbed 
nozzle  would  have  resulted  in  much  wasted  computer  time. 

6.  6  Results  of  Off-Design  Cases 

In  neither  case  was  any  appreciable  difference  seen  in  the 
results  in  the  inviscid  core  near  the  nozzle  exit.  Plots  of  Mach 
number  vs.  radial  distance  and  Pitot  pressure  vs  radial  distance 
for  the  first  case  are  presented  in  Figs.  122  and  123.  Plots  of 
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Mach  number  vs.  radial  distance  and  Pitot  pressure  vs  radial 
distance  for  the  second  case  are  presented  in  Figs.  124  and  125. 
The  largest  change  occurred  in  the  near-wall  boundary  layer  as 
expected,  since  the  remainder  of  the  flow  is  supersonic/hypersonic 
and  shouldn't  be  affected  by  a  downstream  pressure  change  except 
through  the  impression  of  the  boundary  layer.  It  is  believed  that 
the  changes  in  the  inviscid  uniform  core  are  in  large  part  due  to 
the  different  ways  in  which  the  back  pressure  is  specified  in  the 
two  solvers. 

6.  7  Conclusions 

Sects.  6. 1-6. 3  demonstrated  the  fact  that,  given  a  somewhat 
favorable  starting  point  for  the  optimization  scheme,  the  automatic 
code  would  seek  out  the  "global  best”  design  of  a  hypersonic  nozzle 
synonymous  with  the  response  surface  global  minimum  of  the  nozzle's 
computed  composite  error.  Through  three  successful  attempts  to 
locate  the  global  minimum,  the  optimization  code  showed  indications 
of  robustness.  An  indication  of  the  accuracy  associated  with  the 
optimization  code  can  be  inferred  from  the  comparative  results  of 
Ch.  4  for  the  Mach  12  nozzle  (Sects.  4. 2. 3.1  and  4. 3. 3.1). 

Sect.  6.4  demonstrates  that  given  a  relatively  unfavorable 
starting  point,  the  optimization  code  may  seek  out  the  best  nozzle 
design  associated  with  the  local  minimum,  though  this  might  not  be 
the  global  minimum.  Nonetheless,  the  optimization  code  still 
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locates  the  best  nozzle  design  in  the  area  in  which  it  begins  its 
search. 

The  primary  objective  has  thus  been  satisfied  in  that  an 
accurate  and  robust  optimization  code  has  been  developed  which, 
with  simple  response  surface  methods  and  state-of-the-art  NS  and 
PNS  codes,  was  used  to  redesign  the  testbed  hypersonic  nozzle  with 
fixed  values  of  nozzle  length,  throat  radius  and  exit  radius. 

There  are  however  two  caveats  to  the  statement  that  the 
optimization  code  is  accurate  and  robust.  First,  the  designer  must 
have  a  fairly  good  idea  of  where  to  begin  the  optimization  search 
for  the  "global  best"  nozzle,  lest  some  local  best  nozzle  be 
obtained.  Second,  and  of  more  importance,  the  nozzle  flow  is  not 
as  disturbance-free  as  a  designer  might  ultimately  prefer.  It  is 
strongly  suspected  after  discussions  with  Korte  that  this  is  due  to 
the  use  of  only  two  design  parameters  [53] .  Had  there  been  more 
time  allowed,  research  into  alternate  or  additional  design 
parameters  such  as  nozzle  length  or  exit  radius  would  have  been 
conducted.  Such  was  not  the  case  and  this  area  remains  open  for 
further  research. 

Regardless  of  the  caveats,  the  "global  best"  hypersonic  nozzle 
obtained  with  the  two  design  parameters  is  probably  only  within  a 
few  percent  of  the  perfect  "global  best"  nozzle,  that  associated 
with  the  variation  of  all  of  the  design  parameters  necessary  to 
achieve  a  perfectly  designed  nozzle  which  will  provide  a'  completely 
disturbance-free,  uniform,  inviscid  core  at  the  design  Mach  niomber 
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at  the  nozzle  exit.  The  current  optimization  code's  "global  best" 
nozzle  thus  represents  an  extremely  good  engineering  solution  to 
the  hypersonic  nozzle  design  problem,  and  while  not  delivering  the 
perfect  nozzle  flow,  provides  at  the  very  least  a  good  baseline 
nozzle  which  can  be  adjusted  with  further  research  to  provide  the 
perfect  nozzle  flow.  Note  however  that  the  optimization  code's 
best  nozzle  is  obtained  in  a  relatively  small  amount  of  time  due  to 
the  fact  that  only  two  design  parameters  are  used.  Thus,  in  the 
current  optimization  code  trade-offs  may  exist  between  solution 
accuracy  and  time.  The  addition  or  substitution  of  one  or  more 
design  parameters  in  further  research  may  minimize  or  eliminate  any 
such  trade-offs. 

From  Sect.  6.5,  it  is  concluded  that  the  use  of  an  NECR  code 
in  the  optimization  of  the  testbed  nozzle  would  have  likely 
amounted  to  much  wasted  computer  time.  This  is  not  to  suggest  that 
in  the  design  of  a  hypersonic  nozzle  one  may  completely  disregard 
nonequilibrium  and  chemical  reaction  effects.  In  point  of  fact,  it 
is  very  likely  that  these  effects  should  be  accounted  for  in  at 
least  the  last  optimization  iteration.  Though  not  demonstrated 
here,  a  PNS  version  of  the  NECR  code  was  developed  and  validated 
during  the  course  of  this  research  and  is  currently  ready  to  be 
assimilated  along  with  its  NS  counterpart  into  an  NECR  hypersonic 
nozzle  optimization  code.  This  brings  up  the  point  that  the 
optimization  code  as  designed  could  be  run  with  almost  any  NS  and 
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PNS  codes  if  one  paid  careful  attention  to  the  variables  passed 
between  the  parts  of  the  optimization  code  and  the  flow  codes. 

From  Sect.  6.6  it  is  concluded  that,  though  the  capability  to 
test  for  off-design  conditions  is  likely  important  to  the  design 
engineer,  unless  drastic  changes  are  made  in  the  back  pressure  of 
the  nozzle  relatively  insignificant  changes  will  occur  in  the 
nozzle  flow.  The  optimized  nozzle  was  nevertheless  tested  for  the 
effect  of  changing  back  pressure  and  the  results  analyzed, 
satisfying  the  final  objective  of  this  research. 

In  summary,  all  of  the  objectives  for  this  research  were 
satisfied  and  an  excellent  engineering  design  tool  has  been 
created. 
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Appendix  A;  "WL  Mach  12  Nozzle  Wall  Data  (x,y)  (m) 


Nozzle  Wall 


-0,74295E-01 

0.32916E-01 

-0.10160E-01 

0.14526E-01 

-0.73660E-01 

0.32736E-01 

-0.76200E-02 

0.13797E-01 

-0.71120E-01 

0.32007E-01 

-0.50800E-02 

0.13068E-01 

-0.68580E-01 

0.31278E-01 

-0.25400E-02 

0.12342E-01 

-0.66040E-01 

0.30549E-01 

O.OOOOOE+00 

0.11613E-01 

-0.63500E-01 

0.29822E-01 

0.12700E-02 

0.11250E-01 

-0.60960E-01 

0.29093E-01 

0.25400E-02 

0.10884E-01 

-0.58420E-01 

0.28364E-01 

0.38100E-02 

0.10521E-01 

-0,55880E-01 

0.27635E-01 

0.50800E-02 

0.10155E-01 

-0.53340E-01 

0.26909E-01 

0.63500E-02 

0.97917E-02 

-0.50800E-01 

0.26180E-01 

0.76200E-02 

0.94285E-02 

-0.48260E-01 

0.25451E-01 

0.88900E-02 

0.90627E-02 

-0.45720E-01 

0.24722E-01 

0.10160E-01 

0.86995E-02 

-0.43180E-01 

0.23995E-01 

0.11430E-01 

0.83363E-02 

-0.40640E-01 

0.23266E-01 

0.12700E-01 

0.79705E-02 

-0.38100E-01 

0.22537E-01 

0.13970E-01 

0.76073E-02 

-0.35560E-01 

0.21808E-01 

0.15240E-01 

0.72415E-02 

-0.33020E-01 

0.21082E-01 

0.16510E-01 

0.68783E-02 

-0.30480E-01 

0.20353E-01 

0.17780E-01 

0.65151E-02 

-0.27940E-01 

0.19624E-01 

0.19050E-01 

0.61493E-02 

-0.25400E-01 

0.18895E-01 

0.20320E-01 

0.57861E-02 

-0.22860E-01 

0.18169E-01 

0.21590E-01 

0.54229E-02 

-0.20320E-01 

0.17440E-01 

0.22428E-01 

0.51816E-02 

-0.17780E-01 

0.16711E-01 

0.22860E-01 

0.51816E-02 

-0.15240E-01 

0.15982E-01 

0.23368E-01 

0.51816E-02 

-0.12700E-01 

0.15255E-01 

0.23876E-01 

0.51816E-02 

145 


0.24384E-01 
0.24892E-01 
0.25400E-01 
0.25908E-01 
0.26416E-01 
0.26924E-01 
0.27432E-01 
0.27940E-01 
0.28448E-01 
0.28956E-01 
0.29464E-01 
0.29972E-01 
0.30480E-01 
0.30988E-01 
0.31496E-01 
0.32004E-01 
0.32512E-01 
0.32791E-01 
0.33020E-01 
0.35560E-01 
0.38100E-01 
0.40640E-01 
0.43180E-01 
0,45720E-01 
0.48260E-01 
0.50800E-01 
0.53340E-01 
0.55880E-01 
0.58420E-01 
0. 60960E-01 
0.63500E-01 
0. 66040E-01 
0.68580E-01 
0.71120E-01 
0.73660E-01 
0.76200E-01 
0.78740E-01 
0.81280E-01 
0.83820E-01 
0. 86360E-01 
0.88900E-01 
0.91440E-01 
0.93980E-01 
0.96520E-01 
0.99060E-01 
0.10160E+00 
0.10414E+00 
0.10668E+00 
0.10922E+00 


0.51816E-02 
0.51816E-02 
0.51816E-02 
0.51816E-02 
0.51816E-02 
0.51816E-02 
0.51816E-02 
0.51816E-02 
0.51816E-02 
0.51816E-02 
0.51816E-02 
0.51816E-02 
0.51816E-02 
0.51816E-02 
0.51816E-02 
0.51816E-02 
0.51816E-02 
0.51816E-02 
0.52172E-02 
0.56185E-02 
0. 60198E-02 
0. 64211E-02 
0.68199E-02 
0.72212E-02 
0.76225E-02 
0.80239E-02 
0.84252E-02 
0.88240E-02 
0.92253E-02 
0.96266E-02 
0.10028E-01 
0.10427E-01 
0.10828E-01 
0.11229E-01 
0.11631E-01 
0.12032E-01 
0.12431E-01 
0.12832E-01 
0.13233E-01 
0.13635E-01 
0.14036E-01 
0.14435E-01 
0.14836E-01 
0.15237E-01 
0.15639E-01 
0.16038E-01 
0.16439E-01 
0.16840E-01 
0.17242E-01 


0.11176E+00 

0.11430E+00 

0.11684E+00 

0.11938E+00 

0,12192E+00 

0.12446E+00 

0.12700E+00 

0.12954E+00 

0.13208E+00 

0.13462E+00 

0.13716E+00 

0.13970E+00 

0.14224E+00 

0,14478E+00 

0.14732E+00 

0.14986E+00 

0.15240E+00 

0.15494E+00 

0.15748E+00 

0.16002E+00 

0.16256E+00 

0.16510E+00 

0.16764E+00 

0.17018E+00 

0.17272E+00 

0.17526E+00 

0.17780E+00 

0.18034E+00 

0.18288E+00 

0.18542E+00 

0.18796E+00 

0.19050E+00 

0.19304E+00 

0.19558E+00 

0.19812E+00 

0.20066E+00 

0.20320E+00 

0.21590E+00 

0.22860E+00 

0.24130E+00 

0.25400E+00 

0.26670E+00 

0.27940E+00 

0.29210E+00 

0.30480E+00 

0.31750E+00 

0.33020E+00 

0.34290E+00 

0.35560E+00 


0.17643E-01 

0.18042E-01 

0.18443E-01 

0.18844E-01 

0.19246E-01 

0.19647E-01 

0.20046E-01 

0.20447E-01 

0.20848E-01 

0.21250E-01 

0.21651E-01 

0.22050E-01 

0.22451E-01 

0.22852E-01 

0.23254E-01 

0.23652E-01 

0.24054E-01 

0.24455E-01 

0.24856E-01 

0.25258E-01 

0.25657E-01 

0.26058E-01 

0.26459E-01 

0.26860E-01 

0.27262E-01 

0.27661E-01 

0.28062E-01 

0.28463E-01 

0.28865E-01 

0.29263E-01 

0.29665E-01 

0.30066E-01 

0.30467E-01 

0.30869E-01 

0.31267E-01 

0.31669E-01 

0.32070E-01 

0.34066E-01 

0.36048E-01 

0.38014E-01 

0.39967E-01 

0.41902E-01 

0.43825E-01 

0.45733E-01 

0.47625E-01 

0.49505E-01 

0.51369E-01 

0.53218E-01 

0.55055E-01 


146 


0.36830E+00 
0.38100E+00 
0.39370E+00 
0.40640E+00 
0.41910E+00 
0.43180E+00 
0.44450E+00 
0.45720E+00 
0.46990E+00 
0.48260E+00 
0.49530E+00 
0.50800E+00 
0.52070E+00 
0.53340E+00 
0.54610E+00 
0.55880E+00 
0.57150E+00 
0.58420E+00 
0.59690E+00 
0. 60960E+00 
0.62230E+00 
0. 63500E+00 
0.64770E+00 
0.66040E+00 
0.67310E+00 
0. 68580E+00 
0.69850E+00 
0.71120E+00 
0.72390E+00 
0.73660E+00 
0.74930E+00 
0.76200E+00 
0.77470E+00 
0.78740E+00 
0.80010E+00 
0.81280E+00 
0.82550E+00 
0.83820E+00 
0.85090E+00 
0.86360E+00 
0.87630E+00 
0.88900E+00 
0.90170E+00 
0.91440E+00 
0.92710E+00 
0.93980E+00 
0.95250E+00 
0.96520E+00 
0.97790E+00 


0.56876E-01 
0.58684E-01 
0.60480E-01 
0.62260E-01 
0. 64028E-01 
0.65781E-01 
0. 67521E-01 
0. 69248E-01 
0.70963E-01 
0.72664E-01 
0.74353E-01 
0.76027E-01 
0.77691E-01 
0.79339E-01 
0.80978E-01 
0.82603E-01 
0.84216E-01 
0.85816E-01 
0.87404E-01 
0.88981E-01 
0.90546E-01 
0.92098E-01 
0.93640E-01 
0.95169E-01 
0.96685E-01 
0.98191E-01 
0.99687E-01 
0.10117E+00 
0.10264E+00 
0.10410E+00 
0.10555E+00 
0.10699E+00 
0.10842E+00 
0.10984E+00 
0.11125E+00 
0.11264E+00 
0.11403E+00 
0.11540E+00 
0.11677E+00 
0.11812E+00 
0.11947E+00 
0.12080E+00 
0.12213E+00 
0.12344E+00 
0.12475E+00 
0.12604E+00 
0.12733E+00 
0.12860E+00 
0.12987E+00 


0.99060E+00 

0.10033E+01 

0.10160E+01 

0.10287E+01 

0.10414E+01 

0.10541E+01 

0.10668E+01 

0.10795E+01 

0.10922E+01 

0.11049E+01 

0.11176E+01 

0.11303E+01 

0.11430E+01 

0.11557E+01 

0.11684E+01 

0.11811E+01 

0.11938E+01 

0.12065E+01 

0.12192E+01 

0.12319E+01 

0.12446E+01 

0.12573E+01 

0.12700E+01 

0.12827E+01 

0.12954E+01 

0.13081E+01 

0.13208E+01 

0.13335E+01 

0.13462E+01 

0.13589E+01 

0.13716E+01 

0.13843E+01 

0.13970E+01 

0.14097E+01 

0.14224E+01 

0.14351E+01 

0.14478E+01 

•0.14605E+01 

0.14732E+01 

0.14859E+01 

0.14986E+01 

0.15113E+01 

0.15240E+01 

0.15367E+01 

0.15494E+01 

0.15621E+01 

0.15748E+01 

0.15875E+01 

0.16002E+01 


0.13113E+00 
0.13237E+00 
0.13361E+00 
0.13484E+00 
0.13606E+00 
0.13727E+00 
0.13848E+00 
0.13967E+00 
0.14085E+00 
0.14203E+00 
0.14320E+00 
0.14435E+00 
0.14550E+00 
0.14664E+00 
0.14778E+00 
0.14890E+00 
0.15002E+00 
0.15113E+00 
0.15223E+00 
0.15332E+00 
0.15440E+00 
0.15548E+00 
0.15655E+00 
0.15761E+00 
0.15866E+00 
0.15971E+00 
0.16075E+00 
0.16178E+00 
0.16280E+00 
0.16382E+00 
0.16483E+00 
0.16583E+00 
0.16682E+00 
0. 16781E+00 
0.16879E+00 
0.16976E+00 
0.17073E+00 
0.17169E+00 
0.17264E+00 
0.17358E+00 
0.17452E+00 
0.17546E+00 
0.17638E+00 
0.17730E+00 
0.17821E+00 
0.17912E+00 
0.18002E+00 
0.18092E+00 
0.18181E+00 
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0.16129E+01 

0.16256E+01 

0.16383E+01 

0.16510E+01 

0.16637E+01 

0.16764E+01 

0.16891E+01 

0.17018E+01 

0.17145E+01 

0.17272E+01 

0.17399E+01 

0.17526E+01 

0.17653E+01 

0.17780E+01 

0.17907E+01 

0.18034E+01 

0.18161E+01 

0.18288E+01 

0.18415E+01 

0.18542E+01 

0.18669E+01 

0.18796E+01 

0.18923E+01 

0.19050E+01 

0.19177E+01 

0.19304E+01 

0.19431E+01 

0.19558E+01 

0.19685E+01 

0.19812E+01 

0.19939E+01 

0.20066E+01 

0.20193E+01 

0.20320E+01 

0.20447E+01 

0.20574E+01 

0.20701E+01 

0,20828E+01 

0.20955E+01 

0.21082E+01 

0.21209E+01 

0.21336E+01 

0.21463E+01 

0.21590E+01 

0.21717E+01 

0.21844E+01 

0.21971E+01 

0.22098E+01 

0.22225E+01 


0.18269E+00 

0.18356E+00 

0.18443E+00 

0.18530E+00 

0.18615E+00 

0.18700E+00 

0.18785E+00 

0.18869E+00 

0.18952E+00 

0.19036E+00 

0.19209E+00 

0.19200E+00 

0.19281E+00 

0.19361E+00 

0.19442E+00 

0.19521E+00 

0.19600E+00 

0.19679E+00 

0.19757E+00 

0.19834E+00 

0.19911E+00 

0.19988E+00 

0.20064E+00 

0.20139E+00 

0.20214E+00 

0.20289E+00 

0.20363E+00 

0.20436E+00 

0.20509E+00 

0.20582E+00 

0,20654E+00 

0.20726E+00 

0.20797E+00 

0.20870E+00 

0.20938E+00 

0.21008E+00 

0.21077E+00 

0.21146E+00 

0.21215E+00 

0.21283E+00 

0.21351E+00 

0.21418E+00 

0.21485E+00 

0,21551E+00 

0.21617E+00 

0.21683E+00 

0.21748E+00 

0.21813E+00 

0.21877E+00 


0.22352E+01 

0.22479E+01 

0.22606E+01 

0.22733E+01 

0.22860E+01 

0.22987E+01 

0.23114E+01 

0.23241E+01 

0.23368E+01 

0.23495E+01 

0.23622E+01 

0.23749E+01 

0.23876E+01 

0.24003E+01 

0.24130E+01 

0.24257E+01 

0.24384E+01 

0.24511E+01 

0.24638E+01 

0.24765E+01 

0.24892E+01 

0.25019E+01 

0.25146E+01 

0.25273E+01 

0.25400E+01 

0.25527E+01 

0.25654E+01 

0.25781E+01 

0.25908E+01 

0.26035E+01 

0.26162E+01 

0.26289E+01 

0.26416E+01 

0.26543E+01 

0.26670E+01 

0.26797E+01 

0.26924E+01 

0.27051E+01 

0.27178E+01 

0.27305E+01 

0.27432E+01 

0.27559E+01 

0.27686E+01 

0.27813E+01 

0.27940E+01 

0.28067E+01 

0.28194E+01 

0.28321E+01 

0.28448E+01 


0.21941E+00 

0.22005E+00 

0.22068E+00 

0.22131E+00 

0.22194E+00 

0.22256E+00 

0.22317E+00 

0.22379E+00 

0.22440E+00 

0.22501E+00 

0.22561E+00 

0.22621E+00 

0.22680E+00 

0.22740E+00 

0.22799E+00 

0.22857E+00 

0.22915E+00 

0.22973E+00 

0.23031E+00 

0.23088E+00 

0.23145E+00 

0.23201E+00 

0.23258E+00 

0.23314E+00 

0.23369E+00 

0.23425E+00 

0.23480E+00 

0.23534E+00 

0.23589E+00 

0.23643E+00 

0.23697E+00 

0.23750E+00 

0.23803E+00 

0.23856E+00 

0.23909E+00 

0.23961E+00 

0.24013E+00 

0.24065E+00 

0.24117E+00 

0.24168E+00 

0.24219E+00 

0.24269E+00 

0.24320E+00 

0.24370E+00 

0.24420E+00 

0.24470E+00 

0.24519E+00 

0.24568E+00 

0.24617E+00 
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0.28575E+01 

0.28702E+01 

0.28829E+01 

0.28956E+01 

0.29083E+01 

0.29210E+01 

0.29337E+01 

0.29464E+01 

0.29591E+01 


0.24665E+00 

0.24714E+00 

0.24762E+00 

0.24810E+00 

0.24857E+00 

0.24905E+00 

0.24952E+00 

0.24999E+00 

0.25045E+00 


0.29718E+01 

0.29845E+01 

0.29972E+01 

0.30099E+01 

0.30226E+01 

0.30353E+01 

0.30480E+01 


0.25092E+00 

0.25138E+00 

0.25184E+00 

0.25230E+00 

0.25275E+00 

0.25321E+00 

0.25365E+00 
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Appendix  B:  Search  Path  Tabulated  Data 


B.l  Tabulated  Path  Data  for  Automatic  Search  from  Sect.  6.1 

©attach  ©eait  Er  Erl  Ert 

(Begin  first-order  search) 

.1276873E+02  .2291832E+01  .7081753E-01  .9296892E-01  .1637864E+00 

.1776873E+02  .3791832E+01  .1177701E+00  .1003611E+00  .2181312E+00 

.1776873E+02  .7918313E+00  .2453667E-01  .1020390E+00  .1265757E+00 

.7768725E+01  .7918313E+00  .1795622E+00  .1560986E+00  .3356608E+00 

.7768725E+01  .3791832E+01  .5881502E-01  .1098336E+00  .1686486E+00 

(Begin  sampling  in  direction  of  steepest  descent) 

.1436688E+02  .2518504E+01  .2513864E-01  .1702819E+00  .1954205E+00 

.1596504E+02  .2745177E+01  .1965192E-01  .1108700E+00  .1305219E+00 

.1756320E+02  .2971849E+01  .3326964E-01  .9911651E-01  .1323861E+00 

(Inflection  point  detected,  calculate  and  run  case) 
.1345775E+02  .2389559E+01  .5058493E-01  .1578817E+00  .2084667E+00 
(Neighboring  point,  not  gne„,  is  minimum;  begin  second  iteration, 
first-order  search) 

.1596504E+02  .2745177E+01  .1965192E-01  .1108700E+00  .1305219E+00 

.1763174E+02  .3245177E+01  .2584804E-01  .1883572E+00  .2142052E+00 

.1763174E+02  .2245177E+01  .1092248E+00  .1031036E+00  .2123284E+00 

.1429834E+02  .2245177E+01  .1568505E-01  .9243594E-01  .1081210E+00 

.1429834E+02  .3245177E+01  .1102343E+00  .1001207E+00  .2103550E+00 

(Begin  sampling  in  direction  of  steepest  descent) 

.1554069E+02  .2622524E+01  .1907666E-01  .1085871E+00  .1276638E+00 

.1511633E+02  .2499870E+01  .1804721E-01  .1062700E+00  .1243172E+00 

.1469198E+02  .2377217E+01  .1649691E-01  .1105367E+00  .1270336E+00 

(Inflection  point  detected,  calculate  g^ew  run  case) 
.1509428E+02  .2493496E+01  .1801868E-01  .1047688E+00  .1227875E+00 
(gne„  is  new  minimum,  begin  third  iteration,  first-order  search) 
.1509428E+02  .2493496E+01  .1801868E-01  .1047688E+00  .1227875E+00 

.1564984E+02  .2660166E+01  .1998629E-01  .1060464E+00  .1260327E+00 

.1564984E+02  .2326826E+01  .1775123E-01  .9388037E-01  .1116316E+00 

.1453872E+02  .2326826E+01  .1612252E-01  .1096304E+00  .1257529E+00 

.1453872E+02  .2660166E+01  .3686050E-01  .1916600E+00  .2285205E+00 

(Begin  sampling  in  direction  of  steepest  descent) 

.1523284E+02  .2451729E+01  .2239647E-01  .9700075E-01  .1193972E+00 

.1537139E+02  .2409962E+01  .2656149E-01  .9102097E-01  .1175825E+00 

.1550995E+02  .2368195E+01  .2174513E-0I  .8891058E-01  .1106557E+00 

.1564851E+02  .2326428E+01  .1775573E-01  .9413861E-01  .1118943E+00 

(Inflection  point  detected,  calculate  gnew  run  case) 
.1555821E+02  .2353647E+01  .1881112E-01  .8705963E-01  .1058707E+00 
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ig^sw  is  iisw  minimum,  begin  fourth  iteration,  first-order  search) 
.1555821E+02  .2353647E+01  .1881112E-01  .8705963E-01  .1058707E+00 
.1574340E+02  .2453647E+01  .2031832E-01  .8986763E-01  .1101859E+00 
.1574340E+02  .2253647E+01  .3034221E-01  .8751059E-01  .1178528E+00 
.1537302E+02  .2253647E+01  .1741287E-01  .7207244E-01  .8948531E-01 
.1537302E+02  .2453647E+01  .2606265E-01  .9388064E-01  .1199433E+00 
(Begin  sampling  in  direction  of  steepest  descent) 

.1551680E+02  .2326262E+01  .1816967E-01  .8889402E-01  .1070637E+00 
.1547539E+02  .2298876E+01  .1777013E-01  .8367001E-01  .1014401E+00 
.1543398E+02  .2271491E+01  .1752177E-01  .7724496E-01  .9476673E-01 
.1539256E+02  .2244105E+01  .1709301E-01  .7169487E-01  .8878788E-01 
.1534701E+02  .2213982E+01  .1654894E-01  .6469413E-01  .8124307E-01 
.1529690E+02  .2180845E+01  .1662438E-01  .5798683E-01  .7461120E-01 
.1524179E+02  .2144395E+01  .1541290E-01  .5056105E-01  .6597395E-01 
.1518115E+02  .2104300E+01  .1595800E-01  .4277686E-01  .5873486E-01 
.1511446E+02  .2060196E+01  .1565748E-01  .3553693E-01  .5119441E-01 
.1504110E+02  .2011681E+01  .1593548E-01  .2789349E-01  .4382898E-01 
.1496040E+02  .1958314E+01  .1597533E-01  .2160643E-01  .3758176E-01 
.1487163E+02  .1899611E+01  .1635641E-01  .1584001E-01  .3219641E-01 
.1477398E+02  .1835038E+01  .1622345E-01  .1383636E-01  .3005981E-01 
.1466657E+02  .1764007E+01  .1789275E-01  .1607675E-01  .3396950E-01 
(Inflection  point  detected,  calculate  g^gw  and  run  case) 
.1478431E+02  .1841869E+01  .1592740E-01  .1392990E-01  .2985730E-01 
(gng„  is  new  minimum,  begin  fifth  iteration,  first-order  search) 
.1478431E+02  .1841869E+01  .1592740E-01  .1392990E-01  .2985730E-01 
.1488431E+02  .1941869E+01  .1697945E-01  .1845712E-01  .3543657E-01 
.1488431E+02  .1741869E+01  .1731621E-01  .1974095E-01  .3705717E-01 
.1468431E+02  .1741869E+01  .1608365E-01  .1810649E-01  .3419014E-01 
.1468431E+02  .1941869E+01  .2370189E-01  .1644068E-01  .4014256E-01 
(Begin  sampling  in  direction  of  steepest  descent) 

.1479813E+02  .1809325E+01  .1503192E-01  .1442836E-01  .2946028E-01 
.1481194E+02  .1776781E+01  .1498256E-01  .1623601E-01  .3121857E-01 
.1482576E+02  .1744236E+01  .1554608E-01  .1886545E-01  .3441153E-01 
(Inflection  point  detected,  calculate  g^ew  and  run  case) 
.1478810E+02  .1832938E+01  .1581149E-01  .1393989E-01  .2975138E-01 
(New  minimum  at  neighboring  point,  not  gne„,  minimum  d6  reached) 
.1479813E+02  .1809325E+01  .1503192E-01  .1442836E-01  .2946028E-01 
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B.2  Tabulated  Path  Data  for  Automatic  Search/Efficiency  Check  from 
Sect.  6.2 

©attach  ©exit  Er  Erl  Ert 

(Begin  first-order  search) 

.1276873E+02  .2291832E+01  .7081753E-01  .9296892E-01  .1637864E+00 

.1776873E+02  .3791832E+01  .1177701E+00  .1003611E+00  .2181312E+00 

.1776873E+02  .7918313E+00  .2453667E-01  .1020390E+00  .1265757E+00 

.7768725E+01  .7918313E+00  .1795622E+00  .1560986E+00  .3356608E+00 

.7768725E+01  .3791832E+01  .5881502E-01  .1098336E+00  .1686486E+00 

(Begin  sampling  in  direction  of  steepest  descent) 

.1436688E+02  .2518504E+01  .2513864E-01  .1702819E+00  .1954205E+00 

.1596504E+02  .2745177E+01  .1965192E-01  .1108700E+00  .1305219E+00 

.1756320E+02  .2971849E+01  .3326964E-01  .9911651E-01  .1323861E+00 

(Inflection  point  detected,  calculate  gne„  and  run  case) 
.1345775E+02  .2389559E+01  .5058493E-01  .1578817E+00  .2084667E+00 
(Neighboring  point,  not  gne„,  is  minimum;  begin  second  iteration, 
first-order  search) 

.1596504E+02  .2745177E+01  .1965192E-01  .1108700E+00  .1305219E+00 
.1696504E+02  .3745177E+01  .1171590E+00  .9977092E-01  .2169299E+00 
.1696504E+02  .1745177E+01  .4858755E-01  .1949168E-01  .6807923E-01 
.1496504E+02  .1745177E+01  .2016666E-01  .1972702E-01  .3989368E-01 
.1496504E+02  .3745177E+01  .1091028E+00  .9822861E-01  .2073314E+00 
(Begin  sampling  in  direction  of  steepest  descent) 

.1592310E+02  .2394120E+01  .2266946E-01  .9805553E-01  .1207250E+00 
.1588116E+02  .2043062E+01  .5843261E-01  .5190862E-01  .1103412E+00 
.1583923E+02  .1692005E+01  .3954539E-01  .1631110E-01  .5585650E-01 
.1579729E+02  .1340948E+01  .2858088E-01  .9878305E-01  .1273639E+00 
(Inflection  point  detected,  calculate  g^ew  run  case) 
.1584206E+02  .1715721E+01  .4048973E-01  .1600031E-01  .5649004E-01 
(Neighboring  point,  not  g^swr  is  minimum;  begin  third  iteration, 
first-order  search) 

.1583923E+02  .1692005E+01  .3954539E-01  .1631110E-01  .5585650E-01 
.1628367E+02  .2136445E+01  .7259736E-01  .7399603E-01  .1465934E+00 
.1628367E+02  .1247565E+01  .2608620E-01  .1124356E+00  .1385218E+00 
.1539479E+02  .1247565E+01  .2716054E-01  .1917570E+00  .2189176E+00 
.1539479E+02  .2136445E+01  .1966409E-01  .5634728E-01  .7601137E-01 
(Begin  sampling  in  direction  of  steepest  descent) 

.1585064E+02  .1848724E+01  .4862365E-01  .2245980E-01  .7108345E-01 
.1586204E+02  .2005442E+01  .5770983E-01  .4501746E-01  .1027273E+00 
.1587345E+02  .2162161E+01  .5557708E-01  .7452644E-01  .1301035E+00 
(Error  function  not  decreasing  in  direction  of  steepest  descent, 
perform  second-order  search) 

.1583923E+02  .2320538E+01  .2780270E-01  .1027097E+00  .1305124E+00 
.1583923E+02  .1063472E+01  .2760812E-01  .1638776E+00  .1914858E+00 
.1646776E+02  .1692005E+01  .4321268E-01  .1304292E-01  .5625560E-01 
.1521070E+02  .1692005E+01  .3025850E-01  .2397966E-01  .5423816E-01 
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(With  nine-point  stencil  done,  compute  new  minimum  location,  run 
case) 

.1530583E+02  .1954283E+01  .2867062E-01  .2958956E-01  .5826018E-01 
(New  minimum  occurs  on  one  of  nine  stencil  points,  not  calculated 
minimum  point;  begin  fourth  iteration,  first-order  search) 
.1521070E+02  .1692005E+01  .3025850E-01  .2397966E-01  .5423816E-01 

.1550700E+02  .1988305E+01  .4072891E-01  .3911181E-01  .7984071E-01 

.1550700E+02  .1395705E+01  .2926954E-01  .9114892E-01  .1204185E+00 

.1491440E+02  .1395705E+01  .2556720E-01  .1245963E+00  .1501635E+00 

.1491440E+02  .1988305E+01  .1820521E-01  .2244430E-01  .4064951E-01 

(Begin  sampling  in  direction  of  steepest  descent) 

.1520412E+02  .1796556E+01  .2915201E-01  .1974841E-01  .4890042E-01 
.1519754E+02  .1901107E+01  .2635736E-01  .2358737E-01  .4994473E-01 
.1519096E+02  .2005658E+01  .1989707E-01  .3163798E-01  .5153505E-01 
(Inflection  point  detected,  calculate  g^ew  and  run  case) 
.1520798E+02  .1735145E+01  .3005803E-01  .2077856E-01  .5083659E-01 
(Neighboring  point,  not  is  minimum;  begin  fifth  iteration, 

first-order  search) 

.1520412E+02  .1796556E+01  .2915201E-01  .1974841E-01  .4890042E-01 
.1540165E+02  .1994089E+01  .3338242E-01  .3714976E-01  ,7053218E-01 
.1540165E+02  .1599023E+01  .3239598E-01  .3237467E-01  .6477065E-01 
.1500659E+02  .1599023E+01  .2508341E-01  .4107528E-01  .6615869E-01 
.1500659E+02  .1994089E+01  ,1600000E-01  .2581663E-01  .4181662E-01 
(Begin  sampling  in  direction  of  steepest  descent) 

.1514637E+02  .1835824E+01  .2530014E-01  .1954487E-01  .4484501E-01 
.1508861E+02  .1875092E+01  .2073894E-01  .1957144E-01  .4031038E-01 
.1503086E+02  .1914359E+01  .1800037E-01  .1954044E-01  .3754081E-01 
.1497311E+02  .1953627E+01  .1592261E-01  .2178841E-01  .3771103E-01 
(Inflection  point  detected,  calculate  and  run  case) 
.1500533E+02  .1931720E+01  .1648892E-01  .2047684E-01  .3696576E-01 
{gr,ew  is  new  minimum,  begin  sixth  iteration,  first-order  search) 
.1500533E+02  .1931720E+01  .1648892E-01  .2047684E-01  .3696576E-01 
.1513702E+02  .2063410E+01  .1564022E-01  .3569823E-01  .5133846E-01 
.1513702E+02  .1800030E+01  .2606521E-01  .1901789E-01  .4508310E-01 
.1487364E+02  .1800030E+01  .1539135E-01  .1578229E-01  .3117364E-01 
.1487364E+02  .2063410E+01  .2304153E-01  .3118604E-01  .5422757E-01 
(Begin  sampling  in  direction  of  steepest  descent) 

.1498894E+02  .1888139E+01  .1758770E-01  .1726143E-01  .3484912E-01 
.1497256E+02  .1844559E+01  .1702734E-01  .1639507E-01  .3342242E-01 
.1495617E+02  .1800978E+01  .1780645E-01  .1670593E-01  .3451237E-01 
(Inflection  point  detected,  calculate  g^ew  and  run  case) 
.1497146E+02  .1841643E+01  .1727491E-01  .1612557E-01  .3340047E-01 
(gne„  is  new  minimum,  stop  search) 

.1497146E+02  .1841643E+01  .1727491E-01  .1612557E-01  .3340047E-01 
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B.3  Tabulated  Path  Data  for  Automatic  Search/First  Robustness  Test 
from  Sect.  6.3 

©attach  ©exit  Er  Erl  Ert 

(Begin  first-order  search) 

.1500000E+02  .lOOOOOOE+01  .1787525E+00  .4189779E+00  .5977305E+00 

.1750000E+02  .2000000E+01  .2402179E+00  .5937943E-01  .2995973E+00 

.1750000E+02  .OOOOOOOE+00  .2119284E+00  .1134606E+01  .1346535E+01 

.1250000E+02  .OOOOOOOE+00  .5511354E+00  .1884902E+01  .2436037E+01 

.1250000E+02  .2000000E+01  .2253184E+00  .1092264E+00  .3345448E+00 

(Begin  sampling  in  direction  of  steepest  descent) 

.1529728E+02  .1332956E+01  .1730154E+00  .1282010E+00  .3012164E+00 

.1559457E+02  .1665911E+01  .1856333E+00  .1978254E-01  .2054158E+00 

.1589185E+02  .1998867E+01  .2120951E+00  .3823201E-01  .2503271E+00 

(Inflection  point  detected,  calculate  g^e^  and  run  case) 
.1564833E+02  .1726119E+01  .1894446E+00  .1614928E-01  .2055939E+00 
(Neighboring  point,  not  gne„/  is  minimum;  begin  second  iteration, 
first-order  search) 

.1559457E+02  .1665911E+01  .1856333E+00  .1978254E-01  .2054158E+00 
.1609457E+02  .2332581E+01  .2057809E+00  .9362358E-01  .2994045E+00 
.1609457E+02  .9992410E+00  .1683696E+00  .2732370E+00  .4416067E+00 
.1509457E+02  .9992410E+00  .1781558E+00  .4080012E+00  .5861570E+00 
.1509457E+02  .2332581E+01  .1753285E+00  .7093015E-01  .2462587E+00 
(Begin  sampling  in  direction  of  steepest  descent) 

.1562750E+02  .1897489E+01  .1985624E+00  .2482082E-01  .2233833E+00 
.1566043E+02  .2129067E+01  .1914483E+00  .5584736E-01  .2472957E+00 
.1569336E+02  .2360645E+01  .1609997E+00  .7938124E-01  .2403810E+00 
(Error  function  not  decreasing  in  direction  of  steepest  descent, 
perform  second-order  search) 

.1559457E+02  .2608725E+01  .1693448E+00  .1224428E+00  .2917875E+00 
.1559457E+02  .7230972E+00  .1844394E+00  .6543074E+00  .8387468E+00 
.1630168E+02  .1665911E+01  .1931526E+00  .1207620E-01  .2052288E+00 
.1488746E+02  .1665911E+01  .1609893E+00  .2866385E-01  .1896532E+00 
(With  nine-point  stencil  done,  compute  new  minimum  location,  run 
case) 

.1631888E+02  .1821450E+01  .2068326E+00  .2065254E-01  .2274851E+00 
(New  minimum  occurs  on  one  of  nine  stencil  points,  not  calculated 
minimum  point;  begin  third  iteration,  first-order  search) 
.1488746E+02  .1665911E+01  .1609893E+00  .2866385E-01  .1896532E+00 
.1522079E+02  .2110351E+01  .1578574E+00  .3882423E-01  .1966816E+00 
.1522079E+02  .1221471E+01  .1728989E+00  .2016489E+00  .3745477E+00 
.1455413E+02  .1221471E+01  .1628607E+00  .2533742E+00  .4162349E+00 
.1455413E+02  .2110351E+01  .1668066E+00  .4010521E-01  .2069118E+00 
(Begin  sampling  in  direction  of  steepest  descent) 

.1490312E+02  .1821650E+01  .1539040E+00  .1456808E-01  .1684721E+00 
.1491878E+02  .1977390E+01  .1581821E+00  .1942285E-01  .1776049E+00 
.1493445E+02  .2133129E+01  .1703016E+00  .3713999E-01  .2074416E+00 
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(Inflection  point  detected,  calculate  and  run  case) 
.1490404E+02  .1830820E+01  .1542143E+00  .1464542E-01  .1688597E+00 


(Neighboring  point,  not  is  minimum;  begin  fourth  iteration, 

first-order  search) 


.1490312E+02 

.1512534E+02 

.1512534E+02 

.1468090E+02 

.1468090E+02 


.1821650E+01 
.2117950E+01 
.  1525350E+01 
.  1525350E+01 
.2117950E+01 


.1539040E+00 

.1584942E+00 

.1732771E+00 

.1543828E+00 

.1724890E+00 


.1456808E-01 
.3680250E-01 
.5725219E-01 
. 6872501E-01 
.3674111E-01 


.1684721E+00 

.1952967E+00 

.2305293E+00 

.2231078E+00 

.2092301E+00 


(Begin  sampling  in  direction  of  steepest  descent) 

.1491345E+02  .1925499E+01  .1543873E+00  .1639539E-01  .1707827E+00 
.1492377E+02  .2029348E+01  .1629900E+00  .2397354E-01  .1869635E+00 
.1493410E+02  .2133197E+01  .1702498E+00  .3705599E-01  .2073058E+00 


(Error  function  not  decreasing  in  direction  of  steepest  descent, 
perform  second-order  search) 

.1490312E+02  .2240681E+01  .1733245E+00  .5475095E-01  .2280754E+00 
.1490312E+02  .1402619E+01  .1686627E+00  .1147298E+00  .2833926E+00 
.1521739E+02  .1821650E+01  .1753974E+00  .1895371E-01  .1943511E+00 
.1458885E+02  .1821650E+01  .1655162E+00  .1223477E-01  .1777510E+00 


(With  nine-point  stencil  done,  compute  new  minimum  location,  run 
case) 


.1488558E+02  .1876255E+01  .1534002E+00  .1487534E-01 
(New  minimum  occurs  at  calculated  minimum  point; 
iteration,  first-order  search) 

.1488558E+02  .1876255E+01  .1534002E+00  .1487534E-01 
.1503373E+02  .2073789E+01  .1599586E+00  .3069891E-01 
.1503373E+02  .1678722E+01  .1680399E+00  .2583452E-01 
.1473744E+02  ,1678722E+01  .1510056E+00  .2658516E-01 
.1473744E+02  .2073789E+01  .1714617E+00  .2889088E-01 


.1682755E+00 
begin  fifth 

.1682755E+00 

.1906575E+00 

.1938744E+00 

.1775907E+00 

.2003526E+00 


(Begin  sampling  in  direction  of  steepest  descent) 

.1486885E+02  .1810076E+01  .1532397E+00  .1491032E-01  .1681500E+00 

.1485212E+02  .1743896E+01  .1551282E+00  .1854150E-01  .1736697E+00 

.1483539E+02  .1677717E+01  .1571629E+00  .2695767E-01  .1841206E+00 

(Inflection  point  detected,  calculate  g^ew  and  run  case) 
.1487921E+02  .1851063E+01  .1543475E+00  .1448042E-01  .1688279E+00 


(Neighboring  point,  not  g^^^,  is  minimum;  stop  search) 
.1486885E+02  .1810076E+01  .1532397E+00  .1491032E-01  .1681500E+00 
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B.4  Tabulated  Path  Data  for  Automatic  Search/Second  Robustness 
Test  from  Sect.  6.4 

©attach  ©exit  Er  Erl  Ert 

(Begin  first-order  search) 

.lOOOOOOE+02  .3500000E+01  .5476815E-01  .1326148E+00  .1873830E+00 

.llOOOOOE+02  .4500000E+01  .9851435E-01  .1111902E+00  .2097045E+00 

.llOOOOOE+02  .2500000E+01  .1234776E+00  .9904838E-01  .2225260E+00 

.9000000E+01  .2500000E+01  .1190515E+00  .9903943E-01  .2180910E+00 

.9000000E+01  .4500000E+01  .8838512E-01  .1120428E+00  .2004279E+00 

(Begin  sampling  in  direction  of  steepest  descent) 

.9854970E+01  .3822438E+01  .9921345E-01  .1189642E+00  .2181777E+00 

.9709941E+01  .4144876E+01  .7977663E-01  .1177466E+00  .1975233E+00 

.9564911E+01  .4467315E+01  .9363912E-01  .1036490E+00  .1972881E+00 

(Error  function  not  decreasing  in  direction  of  steepest  descent, 
perform  second-order  search) 

.lOOOOOOE+02  .4914214E+01  .1694553E+00  .1594638E+00  .3289191E+00 

.lOOOOOOE+02  .2085786E+01  .1568045E+00  .1227366E+00  .2795411E+00 

.1141421E+02  .3500000E+01  .5570490E-01  .1269973E+00  .1827022E+00 

.8585786E+01  .3500000E+01  .1121790E+00  .1030882E+00  .2152671E+00 

(With  nine-point  stencil  done,  compute  new  minimum  location,  run 
case) 

.9478762E+01  .3455976E+01  .8520553E-01  .9516259E-01  .1803681E+00 
(New  minimum  occurs  at  calculated  minimum  point,  begin  second 
iteration,  first-order  search) 

.9478762E+01  .3455976E+01  .8520553E-01  .9516259E-01  .1803681E+00 

.1014543E+02  .4122646E+01  .9635273E-01  .1133584E+00  .2097112E+00 

.1014543E+02  ,2789306E+01  .6725413E-01  .9792747E-01  .1651816E+00 

.8812092E+01  .2789306E+01  .1401340E+00  .1247406E+00  .2648745E+00 

.8812092E+01  .4122646E+01  .5939336E-01  .1117491E+00  .1711424E+00 

(Begin  sampling  in  direction  of  steepest  descent) 

.9662370E+01  .3603773E+01  .8559036E-01  .9771550E-01  .1833059E+00 

.9845979E+01  .3751571E+01  .9233749E-01  .1087427E+00  .2010801E+00 

.1002959E+02  .3899368E+01  .1073412E+00  .1126093E+00  .2199505E+00 

(Error  function  not  decreasing  in  direction  of  steepest  descent, 
perform  second-order  search) 

.9478762E+01  .4398790E+01  .4761807E-01  .1309440E+00  .1785621E+00 

.9478762E+01  .2513162E+01  .5294066E-01  .6948189E-01  .1224225E+00 

.1042158E+02  .3455976E+01  .2966279E-01  .1328802E+00  .1625430E+00 

.8535948E+01  .3455976E+01  .1147913E+00  .9866041E-01  .2134517E+00 

(With  nine-point  stencil  done,  compute  new  minimum  location,  run 

case) 

.9517540E+01  .3749574E+01  .1074797E+00  .1162802E+00  .2237599E+00 
(New  minimum  occurs  on  one  of  nine  stencil  points,  not  calculated 
minimum  point;  begin  third  iteration,  first-order  search) 

.9478762E+01  .2513162E+01  .5294066E-01  .6948189E-01  .1224225E+00 

.9923202E+01  .2957602E+01  .3594578E-01  .9455798E-01  .1305038E+00 

.9923202E+01  .2068722E+01  .1584257E+00  .1254635E+00  .2838893E+00 
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.9034322E+01  .2068722E+01  .1469209E+00  .8058758E-01  .2275085E+00 

.9034322E+01  .2957602E+01  .1121868E+00  .1440551E+00  .2562419E+00 

(Begin  sampling  in  direction  of  steepest  descent) 

.9555162E+01  .2650472E+01  .4454022E-01  .7150533E-01  .1160455E+00 

.9631562E+01  .2787781E+01  .3931737E-01  .7935443E-01  .1186718E+00 

.9707962E+01  .2925091E+01  .3540313E-01  .9149022E-01  .1268933E+00 

(Inflection  point  detected,  calculate  g^ew  and  run  case) 
.9557502E+01  .2654677E+01  .4437900E-01  .7215990E-01  .1165389E+00 
(Neighboring  point,  not  is  minimum;  begin  fourth  iteration, 

first-order  search) 

.9555162E+01  .2650472E+01  .4454022E-01  .7150533E-01  .1160455E+00 

.9851455E+01  .2946765E+01  .3535811E-01  .9140727E-01  .1267654E+00 

.9851455E+01  .2354178E+01  .1008379E+00  .9185850E-01  .1926964E+00 

.9258869E+01  .2354178E+01  ■.  7955536E-01  .6492129E-01  .1444767E+00 

.9258869E+01  .2946765E+01  .7615678E-01  .1402819E+00  .2164387E+00 

(Begin  sampling  in  direction  of  steepest  descent) 

.9658826E+01  .2635390E+01  .4724712E-01  .8154919E-01  .1287963E+00 
.9762490E+01  .2620308E+01  .5680870E-01  .9303363E-01  .1498423E+00 
.9866154E+01  .2605226E+01  .6760379E-01  .9211896E-01  .1597227E+00 
(Error  function  not  decreasing  in  direction  of  steepest  descent, 
perform  second-order  search) 

.9555162E+01  ,2929822E+01  .4306103E-01  .9608665E-01  .1391477E+00 
.9555162E+01  .2371122E+01  .6928213E-01  .8853766E-01  .1578198E+00 
.9834512E+01  .2650472E+01  .5897811E-01  .8847694E-01  .1474551E+00 
.9275812E+01  .2650472E+01  .6464647E-01  .8534673E-01  .1499932E+00 
(With  nine-point  stencil  done,  compute  new  minimiim  location,  run 
case) 

.9650346E+01  .2738107E+01  .4089450E-01  .7930367E-01  .1201982E+00 
(No  new  minimum,  begin  fifth  iteration,  first-order  search) 
.9555162E+01  .2650472E+01  .4454022E-01  .7150533E-01  .1160455E+00 
.9686849E+01  .2782159E+01  .3886444E-01  .8207822E-01  .1209427E+00 
.9686849E+01  .2518785E+01  .6137357E-01  .8764148E-01  .1490151E+00 
.9423475E+01  .2518785E+01  .5302783E-01  .7527530E-01  .1283031E+00 
.9423475E+01  .2782159E+01  .4859376E-01  .8898575E-01  .1375795E+00 
(Begin  sampling  in  direction  of  steepest  descent) 

.9545297E+01  .2695973E+01  .4306576E-01  .7976903E-01  .1228348E+00 
.9535432E+01  .2741474E+01  .4167120E-01  .8223910E-01  .1239103E+00 
.9525567E+01  .2786975E+01  .4124818E-01  .8468354E-01  .1259317E+00 
(Error  function  not  decreasing  in  direction  of  steepest  descent, 
perform  second-order  search) 

.9555162E+01  .2836710E+01  .4093964E-01  .8668329E-01  .1276229E+00 
.9555162E+01  .2464234E+01  .5910196E-01  .7603489E-01  .1351368E+00 
.9741400E+01  .2650472E+01  .5108793E-01  .8785197E-01  .1389399E+00 
.9368924E+01  .2650472E+01  .5322508E-01  .8570497E-01  .1389300E+00 
(With  nine-point  stencil  done,  compute  new  minimiim  location,  run 
case) 

.9569793E+01  .2692611E+01  .4375831E-01  .7931536E-01  .1230737E+00 
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(No  new  minimum,  begin  sixth  iteration,  first-order 
.9555162E+01  .2650472E+01  .4454022E-01  .7150533E-01 
.9655162E+01  .2750472E+01  .4079748E-01  .8150779E-01 
.9655162E+01  .2550472E+01  .5526874E-01  .8527539E-01 
.9455162E+01  .2550472E+01  .5046561E-01  .7197192E-01 
.9455162E+01  .2750472E+01  .4578677E-01  .8939872E-01 
(Begin  sampling  in  direction  of  steepest  descent) 
.9530786E+01  .2676081E+01  .4371537E-01  .7678889E-01 
.9506411E+01  .2701691E+01  .4396236E-01  .8187671E-01 
.9482035E+01  .2727300E+01  .4468751E-01  .8088665E-01 
(Search  stopped,  no  new  minimum  found) 


search) 
.1160455E+00 
.1223053E+00 
. 1405441E+00 
.1224375E+00 
. 1351855E+00 

.1205043E+00 

.1258391E+00 

.1255742E+00 
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Inviscid/MOC  Inviscid/MOC 


Fig.  2  Departure  of  Actual  Flow  from  MOC/BL  Design  Due  to 
Viscous  Boundary  Layer  Region 


Fig. 


Effect  of  Increased  Edge  Mach  Number  on  Characteristic 
Waves  in  Viscous  Region  of  Boundary  Layer 


Fig.  4  Sivell's  Diagram  for  Wind-Tunnel  Nozzle  Design 


Mach  number 


Fig.  5  Centerline  Mach  Number  Distribution  for  Hypersonic  Nozzle 


Fig.  6  AF  WL  Mach  12  Wind  Tunnel  Nozzle  Profile 
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Fig.  8  Diagram  of  Response  System/Process 
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Fig.  9  Representative  Response  Surface 
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Fig.  10  Representative  Response  Surface  Contour  Plot 


Fig.  11  Steepest  Descent  Method  Five-Point  Sampling  Stencil  in 

Transformed  Space 
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Fig.  16  Size  Reduction  of  Sampling  Stencil  in  Non- 
transformed  Space,  Steepest  Descent  Method 


Fig.  17  Shift  of  Sampling  Stencil  in  Non- transformed  Space  to 
New  Sampling  Stencil  Point  Minimum,  Steepest  Descent  Method 
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Fig.  18  Second-Order  Search  Sampling  Stencil  and  Example 
of  Associated  Response  Functions 


Fig.  19  Transition  to  Next  Sampling  Stencil  in  Non- 
transformed  Space,  Second-Order  Search  Method 
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Fig.  20  Size  Reduction  of  Sampling  Stencil  in  Non- 
trans formed  Space,  Second-Order  Search  Method 


Fig.  21  Shift  of  Sampling  Stencil  in  Non- transformed  Space  to 
New  Sampling  Stencil  Point  Minimum,  Second-Order  Search  Method 
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.  23  Nozzle  Inflow  Boundary  Conditions 
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Fig.  26  Nozzle  Exit  Boundary  Conditions 


Fig.  27  NS  Solver  Provides  Inflow  for  PNS  Solver 
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Fig.  28  PNS  Computational  Grid 


Fig.  29  Grid  Spacing  at  the  Wall 
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n  - 1  n  apace-maiching  direction 

Fig.  30  Determination  of  Wall  Spacing  for  Laminar  Flow 


Nozzle  Wall 


n  - 1  n  space-marching  direction 

Fig.  31  Determination  of  Wall  Spacing  for  Turbulent  Flow 
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Fig.  36  Coarse  53x28  (i,j)  Grid  Used  for  Laminar  Nozzle 

Validation 


Fig.  37  Centerline  Static  Pressure  Plot  for  Laminar  Mach 
12  Nozzle  CFD  Solutions  (53x28) 
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Fig.  38  Centerline  Static  Temperature  Plot  for  Laminar 
Mach  12  Nozzle  CFD  Solutions  (53x28) 


Fig.  39  Centerline  Static  Pressure  Plot  for  Laminar  Mach 
12  Nozzle  CFD  Solutions,  Throat  Region  (53x28) 


188 


1200TI 


1000-; 
^  I 

I  000  T 

T3 

^  600 

e 

L 

• 

?  400- 

K 

200^ 


'  NECR  Cod. .  NECR  Mod. 

I  f  I » r  >  t  >  n  >  yrt t'TrrtTTTT’iniTTrTT  1 1 1  rn  n  i  P'ititti 


-0.10  -0.05  0.00  0.05  0.10  0.15  0.20 

Axlol  Dl«tano»  (m) 


Fig.  40  Centerline  Static  Temperature  Plot  for  Laminar  Mach 
12  Nozzle  CFD  Solutions,  Throat  Region  (53x28) 


P%r  t%c  t  G««  Cod« 


Fig.  41 


Mach  No. 


Profile  of 


Mach  12  Nozzle  Exit, 


Laminar 


Flow,  (53x28) 


189 


1  *-*■■.,*„■*  Per  fee  t  Go#  Cod® 

NECR  Cod®.  P®rf®ot  G«®  Mod* 
tiiif  NECR  Cod®,  hCCR  Mod® 

SAIC  Laminar  CFD  Solution 


I0.  10- 

q: 


0  1 000  2000  3000  4000  5000  6000 

Pr®®®ur®  (  N/ra  2) 


Fig.  42  Pitot  Pressure  of  Mach  12  Nozzle  Exit,  Laminar 

Flow,  (53x28) 
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Fig.  43  Mach  2.0  Flat  Plate  w/Freestream  Data 


Fig.  44  Mach  2.0  Flat  Plate  CFD  Grid 


Fig.  45  CFD  Comparison  of  Viscosity  Ratio  for  Mach  2.0 

Flat  Plate 
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Fig.  46  CFD  Comparison  of  Normalized  U  Component  for 

Mach  2.0  Flat  Plate 


Fig.  47  CFD  Comparison  of  Normalized  Density  for  Mach  2.0 

Flat  Plate 
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Fig.  50  Local  Cf  on  a  Flat  Plate  at  Zero  Incidence  in 

Incompressible  Flow 
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Fig.  51  CFD  Comparison  of  U"  vs.  for  Mach  2.0  Flat  Plate 
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y  (m) 


Fig.  54  Comparison  of  CFD  Results  with  Experimental  Data, 
Mach  8.18  Flat  Plate,  computed  u^  used  in  CFD 


Fig. 


55  Comparison' of  CFD  Results  with  Experimental  Data, 
Mach  8.18  Flat  Plate,  experimental  u^  used  in  CFD 
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Fig.  59  Pitot  Pressure  Comparison  Near  End  of  WL  Mach  12  Nozzle 
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Fig.  62  CFD  Mach  No.  Comparison  Near  End  of  Mach  12 

Nozzle,  B-L  Model  Used 


Fig.  63  CFD  Mach  No.  Comparison  Near  End  of  Mach  12 
Nozzle,  Y-S  Model  Used 
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64  Korte  Mach  2.0  Flat  Plate  w/Freestream  Data 


0.?5  0.50  0.75  1.00 

x(m) 

Fig.  65  Grid  A,  Korte  Mach  2.0  Flat  Plate 
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Fig.  66  Grid  B,  Korte  Mach  2.0  Flat  Plate 
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Fig.  67  Normalized' Temperature  vs.  y,  Korte  Mach  2,0  Flat 
Plate,  GCL  Terms  Not  Included  in  Research  Code 
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Fig.  71  CFD  Comparison  of  Mach  Number  Profile  for  Mach  2.0 

Flat  Plate,  B-L  Model  Used 
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Fig.  72  CFD  Comparison  of  Normalized  U  Component  for  Mach 
2.0  Flat  Plate,  Y-S  Model  Used 


Fig.  73  CFD  Comparison  of  Mach  Number  Profile  for  Mach  2.0 

Flat  Plate,  Y-S  Model  Used 
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Fig.  75  Determination  of  Wall  Spacing  for  WL  Mach  12  Nozzle 
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Fig.  80  Change  in  Pressure  at  PNS  Inflow  Due  to  Acti¬ 
vation  of  Vigneron's  Coefficient  in  Inviscid  Core 
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Fig.  81  Revised  Location  of  PNS  Inflow  Data  in  WL  Mach 

12  Nozzle 
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Fig.  89  Last  Half  of  Nozzle  Separated  into  50  Equal  Regions 
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Fig.  93  Illustration  of  Nozzle  Cone  Angle 


Fig.  94  Illustration  of  Nozzle  Zero  Exit  Angle 
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Fig.  95  Response  Surface  Map  of  Er  Response  Function,  Mach  12 

Nozzle 
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Fig.  96  Response  Surface  Map  of  Erl  Response  Function,  Mach 

12  Nozzle 
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Fig.  97  Response  Surface  Map  of  Ert  Response  Function,  Mach 

12  Nozzle 
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Attachment  Angle  (Degs) 

Fig.  98  Search  Path  Undertaken  with  Second-Order  Search 

Method,  Sect.  5.2. 3.1 


e  LS  Method,  Sect. 
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Fig.  101  Local  Minima  Dete 
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Attachment  Angle  (Deg 

Fig.  104  Search  Path  Undertaken  wit 

Method,  Sect.  5.2.5 


Haoh  No . 

Fig.  105  Exit  Mach  Number  Profile  vs.  Radial  Distance  for 
Optimized  Nozzle  Obtained  Using  Steepest  Descent  Method, 

Sect.  5.2. 5.1 


Fig.  106  C/L  Mach  Number  Profile  vs.  Axial  Distance  for 
Optimized  Nozzle  Obtained  Using  Steepest  Descent 
Method,  Sect.  5. 2. 5.1 
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Fig.  107  Mach  No.  Comparison  Near  End  of  New  and  Old  Nozzle 

Designs  (2.950  m) 


Fig.  108  Mach  No.  Comparison  Along  C/L  of  New  and  Old 

Nozzle  Designs 
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Attachment  Angle  (Degs) 

Fig.  113  Search  Path  Undertaken  During  Automatic  Search, 

Sect.  6.1 
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Fig.  115  Mach  No.  Along  C/L  of  Automatically  Optimized 

Nozzle,  (3.048  m) 


233 


CO  lO  -q-  CO  CM 

(s60a)  0|6uv  ;{X3 


234 


236 


Attachment  Angle  (Degs) 

Fig.  118  Search  Path  Taken  During  Second  Robustness  Test, 

Sect.  6.4 
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Fig.  119  Mach  No.  Comparison  Near  End  of  Newly-Designed 
Nozzle,  Sect.  6.5  (2.950  m) 


Fig.  120  Pitot  Pressure  Comparison  Near  End  of  Newly- 
Designed  Nozzle,  Sect.  6.5  (2.950  m) 
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Fig.  122  Mach  No.  Comparison  Near  End  of  Newly-Designed  Nozzle 
w/Increased  Back  Pressure,  Sect.  6.6  (2.950  m) 
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